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АВ5ТКАСТ 


The aim of this paper is to provide convenient predictor- 
corrector (P-C) methods for obtaining accurate numerical 
Solution at a minimum cost to first order ordinary dirferen 
tial equations (ODE). In pursuing this goal, a unified 
development of the most popular and efficient P-C methods is 
presented, which includes derivation of formulas and analysis 
of error propagation and numerical stability. Each method is 
then coded and programmed using the Fortran language. Com- 
parative analysis of the different P-C methods include both 
theoretical and numerical results. The numerical results 
were obtained by subjecting each method to a wide variety of 
test ODE, using a maximum of two corrector applications and 
a uniform series of step size values. By systematic compari- 
son of the performance of each P-C method the most convenient 


P-C sts in terms of accuracy and minimum cost are established. 
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І. INTRODUCTION 


A linear first order ordinary differential equation (ODE) 
can be used as a mathematical model for a variety of phenom- 
ena, either physical or non-physical. Examples of suelo pi 
nomena include the following: heat flow problems (thermo- 
dynamics}, simple electrical circuits (electrical engineering}, 
force problems (mechanics), rate of bacterial growth (bio- 
logical science), rate of decomposition of radioactive 
material (atomic physics), crystallization rate of a 
chemical compound (chemistry), and rate of population growth 
(statistics). 

Most scientists, engineers or applied mathematicians, 
who have faced the problem of solving an ordinary differential 
equation numerically, are probably aware of the multitude 
of techniques available for such a problem, The abundant 
literature onthe subject of numerical solution of ordinary 
dt tferential eguations is оп (һе опе Папа, а nesult of the 
tremendous variety of actual systems in the physical and 
biological sciences and engineering disciplines that are 
described by ordinary differential equations and, on the 
other hand, a result of the fact that the subject is cur- 
rentlymactive. 

Ое ех Сесе оз засе number of methods, each having 
special advantages, has been a source of confusion as to 


what methods are best for certain classes of problems. It 
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is this observation which particularly motivated Еле тш 
of this paper; thus it is attempted to bring together the 
well-known predictor-corrector methods, which form a class 
of numerical integration methods for solving ordinary dif- 
ferential equations, in a consistent and comprehensive 
framework. 

This delimitation of the study to predictor-corrector 
methods is not without basis. Further research in the field 
сеек а the predictor corrector forms are the most 
efficient among the known integration methods in terms of 
speed and accuracy. Collectively, as a class of integration 
methods, the predictor-corrector sets are the best, but 
maay idually as a predictor-corrector set, the choice tom 
the best method varies depending on the application. This 
paper atempts to compare the overall efficiency of all the 
well-known and most popular predictor-corrector sets by 
using a standard mode of application, the same series of step 
size values, and a set of test ODE's with many unusual and 
interesting features. To introduce the numerical experimen- 
tation of the different predictor-corrector methods, a 
comprehensive analysis of each P-C set starting with its 
derivation, error propagation, and stability is presented 
in detail. 

Numerical experiments are conducted using the IBM 360/67 
digital computer. Тһе tremendous computational capability 
and speed of this computer offered an indispensable tool in 


conducting the experiment on a wide variety of test ODEs. 


PI 





Single precision has seven digit accuracy гап Тош оно 
cision has fourteen digit accuracy available for computation. 
The paper starts with the description of the nature of 
ODE initial value problem. Then the various P-C methods for 
obtaining numerical solutions of the problem are enumerated. 
To lead into the derivation of formulas a brief review of 
backward difference operator and the well-known Newton Back- 
Шара Formula is presented., In turn the different уре ек 
P-C methods are discussed followed by the analysis of error 
propagation and numerical stability. The stability bounds 
of the P-C methods were established through numerical 
experimentation. The next step 1s to subject each P=C 
method to a wide variety ot test ODE. Then systematic com- 
parison of the performance of each P-C method is made.  Con- 
clusions are derived and recommendations are made based on 
the analysis of numerical results obtained. Flow charts and 


computer programs are attached. 


Ша 





Il. NATURE OF -THESPRODPENM 


A. THE NUMERICAL PROBLEM AND NOMENCLATURE 
A linear first-order ordinary differential eguatron 


(ODE) of the form 

Bee a) ВЕТ) 
with the initial condition 

M ML (152) 


where f(x,y) indicates a differentiable function of the 
variables x and y is commonly referred to as an initial 
value problem. | 

A typical elementary differential equations text presents 
several general classes of methods for solving a linear 
first-order ODE. The principal classes of methods are 
(1) variables-separable, or reduction thereto; (2) exact 
equations, or reduction thereto; and (3) solution by infinite 
series. The student is taught to apply the general method 
that appears best for the solution of the particular ODE. 


For example, the linear first-order differential equation 
dy/dx = xy (1-3) 


can easily be solved by the variables-separable method. 


This is accomplished by rewriting the equation in the form 
dy/y = xdx 


and=integrating both sides to obtain 
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1пу = х2/2 е 


where c is an arbitrary constant of integration. пепео пе а 
solution can then be written as 


2 
_ X e 
y = ce 
where су = ес. Тһе general solution of such a linear first- 
order ODE consists of a family curves called the integral 
ox’ /2 


curves. The family of integral curves y = that con- 


T 
stitute the solution of equation (1-3) is illustrated in 


Ре“ 1. 


Figure. l. Solution Curves of y' = Xy. 


lowsweach positive value of cj a particular member of this 
family of curves is determined. 

A particular solution of equation (1-3) can be deter- 
mined if a condition on the solution curve is specified. 
For example, if it is required that the particular solution 


curve pass through the point (0,1), given by the initial 


condition 


y (0) = Yo = Ше 
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22 
then the particular solution y = e^ /? is obtained since 


Si 
However, in real-life problems, many differential equa- 
tions that are encountered cannot be solved by elementary 
classical methods. In such instances, one must resort to 
numerical methods for obtaining one or more particular 
Solutions of the initial value problems A тие о 
techniques are available for solving such a problem numeri- 
cally. This Dane considers only numerical methods using 
predictor-corrector equations. Even with this restriction, 
a large number of methods are in existence. Literature 
ЕБ Сасна this subject; however, reveals that the moct 


popular and efficient methods are the following: 


l ийел Prodictor-Gorrector Method 

2. Milne Predictor-Corrector Method 

3. Nystrom Midpoint Predictor-Euler Corrector Method 
4. Hermite Predictor-Milne Corrector Method 

5. Milne Predictor-Hamming Corrector Method 


6. Adams Predictor-Corrector Methods 


These methods, each having special advantages, have been a 
source of confusion as to what methods are best for certain 
classes of problems. This paper will attempt to present a 
comparative analysis of these various predictor-corrector 
methods. 

In attempting to compare the various predictor-corrector 
methods, such questions as the importance of the number of 


function evaluations, the step size to be used in a numerical 
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calculation, the computing time required Torta шее 
the influence of truncation and round-off error and the 
stability of various predictor=corrector modes toco ае 


tion и111 be considered. 
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III. DERIVATION OF PREDIECIOR-CORREETORZEOUF FRI - 


А. BACKWARD DIFFERENCE OPERATOR 

To lead into the development of certain equations of 
major importance in numerically solving ODE, a linear back- 
ward difference operator V is defined by the following 


equations 
Ои = Y -Y 3 (2-1) 


prc от аз, е 2 


= № ma gi ў 
УТ РЕКУ T A а ТЫ: 


q l! gel 
T E Уһ | | Yn-1 


< 
© 
“2 
1 


ш репета1‚ these formulas start with a given sequence но 


the ЕЕ E the ordinates at equally spaced Xe MAE nh, 


hugs a constant. 


B. NEWTON BACKWARD DIFFERENCE FORMULA 
Using these operators the well known Newton backward 


difference interpolation formula is: 


р 
Мн. IA E (29) 


* уа = т 


Ynta on 





‚ «(ату (оаа) gn 


where a = 
n! Yn? 
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This formula is obtained by fitting a polynomial to the set 
of points (х,у) at equally ѕрасеа {x}. (Actually 

the formula is a polyn mial in a because of the change of 
variable.) Fitting is taken to mean that the polynomial 
passes through the points (X, ,Y4)- Since there are 

mil) points A E this means that if the solution 
function y(x) is a polynomial of degree n or less, the 
formula above will be exact. When y(x) is not a polynomial, 
or is a Perona of degree greater than n, the formula 
will not be exact (except at the points themselves). In 
other words an error will be made because a finite rather 
than an infinite series is used. In general form this error 


will appear as 


т = eh ЕУ 


+1 (Е) is evaluated at 


where C is a function of a and y 
some Ё, X << x. Ty is commonly referred to as the 


Pruncation error. 


C. INTEGRATION FORMULAS FOR ODE 


Integrating both sides or equation (1-1) between (X, 4) 


and СУ а) there results 
Xp +1] 
Уп YA ES Еб,Уах (2-3) 

Zn 

А + 1 

ри y'(x)dx 

Xn 

1 


y *h/ y'(a)da 
i 0 
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Attaching the index sequence of n.= 0,1,2 аға 
tion and noting that Ya is known at Xo? (1-2), it becomes 
apparent that (2-3) can be used recursively to generate 
Y1+>Y2>Y 3»... as long as there is some way to evaluate the 
integral. In other words the solution of the ODE is con- 
verted to evaluating an integral. 
Derivations of many of the equations used to represent 

(2-3) now follows: 

| a) Finite-difference table of backward differences. 
Assuming that (хъ,У.), (х,у)... (х Ур) ате ріуеп, апа 
computing y' = f (X:,Y1) for i - 0,1,...n the following table 


is constructed. 


1 
Хо Yo 
а 
x y! v? ! 
1 1 1? 
уу» 
1 
2 72 Ф 
а 
• V уч 
x ! y? ! 
ІІ 7-1 Уһ 
VS 
Xn Ут 


These values will now be used to continue the solution in 


computing E In general, by retaining gen differences 


ШЕ 





in the Newton Backward Formula (NBF) (2-2), the following 


finite series results. 


-2 ССО ee 
У зо У O E 


: а(о417..0( 54-1) A) yy? (2-4) 


Using this interpolating polynomial to approximate the inte- 
ртапа xn (2-5) [replace y'(x) by Ve yields: (complete 


details of integration are presented in McCalla (ВеҒ. 11) 


4 i | 
DA Sy uU h ET а; У ys with a, 1 
where 
: a(a+l) (ati 1) 2 
a, = J От О (2-5) 
0 1. 


Equation (2-4) is in extrapolation form as can be seen from 


the table of differences, where the end point т 15 


excluded from the interpolating points. 
The error term associated with truncating after the q th 
V 1S 
1 
Э 0 ООСО ории 
а / (ТО Й (Е) de 


[д +2] 


Е ИЕ тре сое Ете теп ора y does not change in sign 


in [0,1], it is possible to write 
бе? [ше] 2 
атты (2-6) 


By actually calculating the а. 1n (2-5) one arrives at 
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B 1 5 2 Хр 
Ya = ЕИ V DPR ИС“ tas У Баја 


Note that by truncating this series at the first difference 


then (for q = 0), 
EP 0-12 
У 1 У, + By, with T, = 1.n^ yl*1 qo . (aap 


which is the well known Euler's formula. Adopting the con- 


vention that 


Та = T(x,h) 


since a is a function of x and h. Also 
T(x,h) = $ nf yl] çg) 


can be written as 
" 2 
MESA) О (а) 


indicating that the truncation error is of order na 
Actually (2-3) can be generalized by rewriting 1t in the 
form 


Jl 


= ! | - 
Ул+1 Yn-r ^ u 6 Ynta ae a 


where г 15 апу positive integer. For example, the case r = O 
15 the one already discussed. Following the same procedure 
as previously presented, the results for r = 1,3,5 are 


Obtalned by substituting (2-4) into (2-8) where 


1 Е 
a. =f alori a s da, f 


1! О =. 


E 


Actually calculating a; for different мајџез ов о ено бе 


results 


Zi 





ч 1.2. 1.3 Ө 
г m. xs Y MS h[2+0V+ z үз 7 V + 30 LES (2-9) 
DET 
Ж” Ж Во о 4 
т SE a ура ре ШЕ 7 V *0V^- 15 У + (2-10) 
Ер l. 
мы _ | f... 5 
T 635 Е gt h[6-12V+15V -9V "+ To ү чору S CUN 


The error term associated with truncation of these formulas 
is more complicated to calculate because the integrand 
changes sign in [-r,1]. For special cases however, as in 
the above for r = 1,3,5, the coefficient of the m dif- 
ference is always zero. Thus it is only necessary to use 
ищет | differences ot y' in order to calculate a result 


mrose truncation error ls of the order of r+2 in h. For 
= ! | = 
У еди | (2-12) 
with 


2 
T(x,h) = 2 y 5 (2) 


Equation (2-12) is known as the Nystrom Midpoint formula. 


Similarly for: 
1935 9735 Y 7 Yu s * 4Һ[у!-Уу) + 5 Уут] (2-13) 
тү ne n n 5 n 


with 
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T(x,h) » li p^ ШЕ 


+ 61 (У = Ша + >. vy! -3 voy! + 11 v^y!] 


У +] 7 7 2 20 п 


n-5 
with 


тб.) = уд ау) 


Generally, formulas derived using the extrapolation form of 
the NBF are referred to as open, explicit or predictor 
equation because Yn+] occurs only on the left hand side of 


the equation. In other words, y сап рек са сп каке 


ntl 
directly from the right-hand side values. 
b) Extending the table of backward differences by one 


line to include Х +] а$ ап interpolating point yields 


x ae 


О О 
! 
VY] 
x У} уу; 
1 1 2 
ү у! 
qu! 
; ү У +] 
vn 
X M : 2 
n n V ЖІ 
уЎп+1 
x 


! 
mel Yn+1 


NBF now in interpolation form yields amfinin ee 


го оше: 
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1 ха ü 1 (e 1) Zu 
Улъа na AA е. Yn«] * 
(a-1) (a) (041) Se (оса NES _ 
+ а! x (2-14) 


Repeating the procedure used in the extrapolation form 


(2-4), analagous results for r = 0,1,3 and 5 are 





т = 0; уд УЗ ше 5 7- 2. 7 = Ar Lv 
УЫ (2-15) 
т = 1; AAA š h[2-2V * š VENUE = = v^ 
Ec (2-16) 
r= 3; a Ya t h[4-8V + v? = š v^ + Ir v^ - 0v? 
·. „Ју +] (2-17) 
r5; y 7 Yy.g * hl6 - 18V « 27v^ - 24v? « 153 g^ 
- eV yl CAs} 


Some interesting and useful results may be obtained by trun- 
cating after a certain number of differences. In the r = 0 


case truncation after the first difference (q-1) yields 
= у + 1 [y! y! Co 
Уп+1 Уһ p [Ул+1 Уһ! 
with 


3 
T(x,h) = “BS y !% (5) 


24 





which is called the modified Euler formula cee 
truncating after the second difference (note лава не enim 


difference is zero) yields 


h 
t ЕС За ТЫ "a ССС 
with 
току) + 289 lS) gy 
2 90 у 
nich is popularly known as Milne's corrector. Formulas 


resulting from the use of NBF interpolation form are referred 
to as closed, implicıt or corrector equations since И 
occurs on both sides of the equation. In other words 


unknown У +] Cannot ре calculated directly since ІІ Ше еселі 


tained within E 


c) In sections a and b the equations were obtained 
directly from manipulating interpolation or extrapolation 
polynomials using NBF. However, another method, called the 
депоа ост Undetermined coefficients, applied to the gen- 
eralized lincar, K-step differential- difference equation 


with constant coefficients 


MA 


K 
y = ЕРУ st has 
n+l ое em) 


will yield all the formulas derived so far and at the same 
time can be used to obtain all other predictor-corrector 
equations by suitable choice of the parameters a; and В; - 
But since the previous formulas derived are sufficient for 


applications and only the Hermite predictor т ets яи 


2 





using the present method, the tedious and Топ рос ЕЕ 
derivations, which can be found in most numerical analysis 
books [Refs. 1,2,3], are left out and only the result mi on 
the Hermite predictor from (Ref. 4] are presented. There 


ise Was found that 


dor -4 a, = Š а, = а, = 0 
Ес 4 B, = 2 Кылы суты s: 0 
yielding 
re E * 5У,-1 XD Ша + 2У1-11 (2-22) 


with 
h^ [4] 
T(x,h) = zu (&) 


d) It is of interest to note that the numerical methods 
developed for solving the initial value problem differ in 
their requirements of the number of starting values. For 
example, Euler's formula (2-7) needs only the initial con- 
dition ACTO ах and h to start the solution. Methods that 
determine Yn+1> when only one point (СУ) and step size 
h are known, are commonly called one-step method. 

On the other hand, the Nystrom Midpoint formula (2-12) 
needs starting values SS and Ce to continue 
the solution. Methods that require step size h and more 
than one point (X, Ya)» (х1-19Ур-1)5:-: in order to compute 
У +] are called multi-step methods: 


It is clear then that a multi-step method = require шя 
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one-step method to provide the necessary starting values 
in its computation. 

Literature research showed that the Runge-Kutta method 
is the best among the available one-step methods. Hence the 
fourth order Runge-Kutta method based on the Taylor's series 
expansion truncated after terms of the fourth order (deri- 
vation of this method is given in Ref. 3) was used as a 
starting method for the multi-step methods considered. This 
has the form 

K1*2K,*2K.*K, 


where 


n” n 
K; à hf (x), + > Yn + = 
= p hE(x, t 2, Yn * = 
Ky = hf (x, ¡A у к.) 


to solve the initial value problem 
SA е 


yielding as many starting values as needed. 
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IV.  PREDICTOR=GORRECTORSMETANOS 


A DEFINITION 

Algorithms that use both an explicit formula and an 
implicit formula are called predictor-corrector methods. 
The solution approximation computed by the explicit formula 
iS denoted c and is called a predictor. This predicí(or 


is then used initially in the right side of an implicit 


C 


formula that computes a corrector MUS The implicit 


C 


За from the 


Formula can be repeatedly applied, using y 
preceding iteration on the right side and computing a new 
M Gini lert.: 

To illustrate the procedure, the Euler P-C algorithm 
eme step method predictor) and the Nystrom P-€ algorichm 
(multi-step method) are used. 


Given the ODE 
у! = f(x,y); with initial solution point (x,>Yp) 


1) Euler P-C method uses the equation (2-7) as predic- 


mer, and (2-19) as corrector 


Been vt is (2-23) 
=y +B [yt _ + y'] (2-24) 
Уп+1 п 2 n+1 n 


Since (2-23) requires only one solution point (x, »Y,) and 
step size h (assumed to be constant for specific application) 
which can be chosen arbitrarily (the choice actually affects 


the stability and error propagation of the method as will 


\ 
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be illustrated later), И. < E can be evaluated. Hence 
225) can be computed yielding ther cult уге (p means 
predicted value). Then the derivative of ун is computed 


yielding the result Уш: 


Е р 
О ST 


Next (2-24) is used to compute a corrected value called 


с 
+ (c means corrected value) 


С h 1 1 | 

= + = + 
Уһ-1 Yn 2 Ua Ху ] · 
At this point there are two alternatives. First, the compu- 
tation can be terminated using и as the true approxima- 
tion value for УН: Then the computation continues to 


determine the next solution point y by repeating the 


11327 
mire procedure. Second, the corrected value, x of 

У +] may not be acceptable in the sense that the equality 
ОТ both sides of (2-24) is not satisfied to as many digits 
as may be desired. In other words, the corrector has not 
converged to the desired accuracy (convergence of the 
corrector formula is an important topic and discussion of 
Ше aspect is deferred to a later section). The next step 


then is to improve the value Ул on the left side of (2-24) 


+1 


by taking the derivative of y calling the result Ур? 
and then calculating an improved value of y This 1s 


ne 


mepeated until convergence occurs сома пало Ба оны ше 
desired. This convergence term, designated as EPS, Ts 
tested against the absolute difference between the last two 


approximations. 
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2) The Nystrom P-C method has (2-12) ас рте ы 


(2-19) as corrector: 


Ve Gawa друг (2-25) 


_ һ _ 
Yn*1 ^ Yn * 2 us y у, | (2-26) 


It is clear that before (2-25) can be evaluated another solu- 
tion point ORE is needed in addition to the given 
solution point. (x ,y, i) of the ODE. Another one step method 
is therefore needed to provide the additional starting 
values. The Runge-Kutta method can be used to accomplish 
this. Once these starting values are known, the same pro- 
седите as in Euler P-C method is followed to compute the 

Next point using equation (2-25) as predictor and (2-26) as 


COrrector. 


EA SIMPLE PREDICTOR-CORRECTOR SET 

The predictor-corrector methods discussed so far do not 
take into account the truncation error incurred in retaining 
th et rereneces of the infinite series of Ес СИРЕ 


application is straightforward, predict then correct, and 


d such it can be called a simple predictor corrector cerr 


C. MODIFIED PREDICTOR-CORRECTOR SET 

The P-C set given by (2-12) and (2-19) has a local атта 
cation error of the same order as the set (2-23), (2-24) 
(T(x,h) = 0(n?)). Hamming [Ref. 4] has used this feature to 
expand the method of calculation. The local truncation 


errors for the predictor (2-12) and corrector (Aim 
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Jr] D 
where 
ÓN Sp EE 
and 
T(x,h). - T - ho ЗЕ) 
С nz lee фа y C 
where 


PA < ох 
n S c n+ | 


Lt follows that the exact value Ve Gs Of yd tx 


given by 
p h^ [3] 
Yxa) T Yn t z Y (6р) 


Э 


from these equations, one obtains 


Са 


5 
с һ [3] u h [5] 
Ма на "ње ЕЙ 


OT 


Са 


5 
ћ 
Уа суы а уа) ту (Е) 


Now assuming that 


151 s ey Е D 
У lp) Уу Созу 6) 


(2229) 


(2 26) 


S 


Gers 


Over the interval of interest (1.е., the third derivative 


does not vary greatly over the interval), (2-29) becomes 


с лы За М5 
Yn«] 7 v er ey, Al) 


Sl 


(AN 





мһеге 


x <. Ë `< x 


n- 1 ТІЛІ 


Even at this point in the development some interesting fea- 
tures can be identified. First one notices that (2-30) can 
be obtained only in case the order of the predictor equation 
is equal to the order of the corrector equation Second: 
there is now a measure of the local truncation error in terms 


and yP which are explicitly calculated in the P-C 


С 
of Уһ ПЕРІ 


+l 


algorithm. This may be seen more explicitly by comparing 


(2-27) with (2-30). 


с ЖА ис аре 


Xu Wt aC) (2-51) 


1 3 
G hay 


3 


| 
jan 


(&)) 


n 


+ | сл 


In*l,p 


Euri; comparing (2-28) with (2-30) 


С 2 VE ines. я 
Улат ~ Ynt1 Я ~> Пере (2-52) 


Since the local truncation error can now be estimated, the 
next question is how to use this information to improve the 
predicted and corrected values. Considering the predicted 


value first, it can be seen from (2-31) that 


= ONES UD 
lirip Б Ол Grae 


Assuming that the local truncation error remains approxi- 


mately constant over two steps, then 


о 





"OR а a: MEN OU В 
In*tl,p 5 У +] Yn+1) 5 Vr E 


The predicted value can then be improved or modified by 

š 4 C я : . 
adding the term [3 Ша = УР) tO у a using only information 
calculated previously: 

The corrector can be modified in essentially the same 


manner since 


E Р 
S 5 У +1 ES 


can be added to y» to improve this value. On this basis 


f 


the overall P-C calculation can be written in the following 


algorithmic steps: 


RE = ! 
Predict: т | Жү 2Һу, 

M s UM ы 
Modify: Пеле pr т с пра Ca) 
Reevaluate ' “ú 
Derivative: 51 f (Xa+? Mp +]) 

= Dom: 
Correct: Che] JE > (m, +1 EL) 
Modify: yt = c E (p АС ) 
` Jt] ns 5 nti n+] 


mem iterate the corrector to convergence. To simplify the 


presentation, the symbols Mi+1? Pei and c have been used. 


n+1 
The use of the моча Ес рессор corrector ea een 

to be very attractive but it must be noted that this method 

mas two limitations. First, it rests on the fact that the 


eriginal P-C equations have equal ова и Еа ико ни vor 


and second, it neglects round-off error. 
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D. HAMMING MODIFIED PREDICTOR-CORRECTOR SET 

Hamming used equations (2-10) and (2-16), which are known 
as the Milne P-C set when both are truncated after the third 
differences, and he revised the corrector to remove the severe 
stability problems, which will be seen later. 

To develop the revised corrector, Hamming started with 


the generalized equation (2-21) and obtained 


= ! ! 
Yn+1 ООД Ы mc Е ne? ч ВВ Yu] T Bo’n 


+ В. Ул 11. 


Using the method of undetermined coefficients [complete deri- 


vatlons can be found in Ref. 4), he obtained 


И OE 
пет (9 да.) 8 1 Т (9 Q4) 
а. = а В = L (9 + 7a.) 
1 | О l: 1 
жо Wi еы шш о. В 
2 8 1 1 24 1 
with 
(-9 + 507) 5 
— ү 


TE) 360 


with a, as a free parameter. Hamming then considered the 
values of а2>92»... 81 as functions of q ТОТ a, = Т, ЭУ 
Ш/9, 0, -1/7, -9/31, апа -6/10. After looking at the result- 
ing coefficients in terms of equal magnitude, number of 
zeros, etc., the case a, = O was chosen as the best value 
since it yields a greater region of real negative stability 


and a wide range of relative stability. It can be verifled 
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р ма ы” C 


that the case a, = 1 yields the Milnes correc ore 


1 
truncated after the third difference. With Q4 7 0 the 


Hamming predictor-corrector set 1s of the Corm 


p n 4 (521 i | Е 
Уһз1 EE 3 n ФА ole ee (2-53) 


3, q = 3, of the NBF extra- 


сл is (2-10) іп Не сазе г 


polation form, with 


15205 [5 
T(x,h) = {8 Һ?у!?!(Е) 
апа 
C | _ 5 ' lay! ра 
vta D wee и (2-34) 
with 
-h? 151 
T(x,h) = -jg Y (8) 


Next, using the procedure discussed in the previous section 
onthe modified predictor-corrector set, Hamming develops 


mrs modified predictor-corrector algorithm as follows: 


eres EX 4 E ' - 
Predict: Pn+1 A 3 h [2y\ Ул-1%2У1-2: (2 55) 
КЕ и 3 me _ 
с ™n+1 > Pn+1 Ga On 
Reevaluate i = £(x га ) 
Derivative: n+1 n+1”° n+1 


= шиш = © t l Zar! 
Correct: Ca+1 О p nh ns и 


€ 9 
/n*1 — n+1 127 Ub cus 





Modify: С ] 


nt] 


Then the corrector may be iterated to convergence as desired. 
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E. P-C SETS CONSIDERED IN THE МОМЕКТСАВНЕ НЕ И Шо 

Having set forth the necessary foundation for the pre- 
dictor-corrector equations, a list of the P-C sets, which 
were analyzed and actually applied to the numerical solution 
of the initial value problem in (1-1) is now summarized. Іп 


each case the operator is evaluated in terms of y/,7» Y» 
1 
Eget. 


EC TI: -Euler P-C set. 


The predictor is equation (2-7) and the corrector is 


(2-24), 
212 
LE EE uM T(x,h) = 5h ye) 
И не u” Е 
ж)” Мәсі ty Iona отра та RN 


P-C-II: Milne P-C set. 
3 of the NBF 


Egquataon (2-10); in the case гт =з, q 
extrapolation form, is used as the predictor, and the cor- 
meetOr 3s (2-16) where r = 1, q = 3, of the NBF interpolation 


form. This yields the equations 


P = 


4 _ | _ 14 „5 [b] 
Your Yous za о У-У ЕСЕ си 
C 2 seo [ 1 +4у! +у! 1: T(x aa Udo 
Yn+1 Yn-1 7 (71-1 n  n+1 у 90 У 


ШӘС-ІТІ: Мувсттопт Ртеайтстог- ЕШ Тек Corrector Set: 
The Predictor is (2-12) and the corrector ea 229 Ән 


yields 


3 
h 3 
УР 1571-15 291: Th) кү 
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C 4 В Г: o 2 ПЕЕ 
Yn+1 "Әл 72 Әні ml ТО, Се 


.Р-С-ТУ: Hermite Predictor-Milne Corrector. 


The predictor used is (2-22) and the corrector is (2 16) 


ylelding 


P gs 
И п Ет 


n-1 


4 
T(x,h) = h^ yl) (gy 


6 
C m h t t 
ОЗЕРЕ PI t un i Ди 
-һ? 
T(x,h) = 2 yg) 


P-C-V: Hamming Modified Predictor-Corrector Set. 
Equation (2-33) is the predictor and (2-34) îs the cor- 


ШЕСТОГ. 


Е 


41 2 
/9+1 E ae МА M zm e 2) 


T(x,h) = те Һ?у!?1(Е) 


c == 
У +] 


coj = 


3h | 
КООР ЛЛ ОЕ сы 


T(x,h) = — > 9131 (5) 


Adams Form (Adams-Bashforth Predictor-Adams- -Moulton Corrector 
t). 


The A-B formulas are (in the case r = 0) equations of 
the NBF extrapolation form and the A-M formulas ате ъп ее 
fac 1-0) equations of the NBF interpolation form. Three 
methods are considered: the second, third, and fourt mordon 


A-B predictor and A-M corrector. 


D 





Р-С-УІ: Second Order ази e Ci 
The predictor and corrector are Obtained sine ue 


q-1 yielding 
DEN h B 
У Ya а у 


T(x,h) = > hyl”! (£) 


C _ h. 
Yn«1 л 2 DES T DM 


13 
аа = е шо 
12 


EO VII: Third Order Adams P=-C Set. 
The case q = 2 in both the extrapolation and interpola- 
tion of the NBF form are used as the predictor and corrector 


equations respectively, resulting in the formulas 


Dine, 2 m и 
о ee о 
9 4 [4 
T(x,h) = g hyt" E) 


€ _ һ | 8 
Yn«l 7 Yn "то [ЗУпњ] t 8Yn T Yn- 


+ 4 
T(x,h) = Bo yl (2) 


Ber VIII: Fourth Order Adams P-€ Set. 
The case q = 3 in both NBF forms yields the predictor- 


corrector equations 


р = h: ! = ! ! S ! 
Yarl Ул za (9597 ° 9911 ФГ 
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T(x,h) 


C 
Yn+1 


T(x,h) 


о 
720 ћу (Е) 


h 


“ y; t 24 Dy tL Е 


From the P-C sets above an interesting fact is observed. 


Namely, P-C-I, 


Реса and P-C-IV used different predictor 


equations but have the same corrector equation while P-C-II 


and P-C-V used the same predictor equation but have dif- 


ferent corrector equations. The significance of this 


observation will be clearly demonstrated later on. 


En 





V. NUMERICAL STABILITY OF PREDICTOR CORRECTOR НЕ 


In the numerical solution of an ODE, a sequence of 
approximations Уһ to the true solution 15552 1s generated. 
Roughly speaking, the stability of a numerical method refers 


to the behavior of the difference or error, у р y x as 
becomes large. In order to begin the discussion, the 
various types of errors which are incurred in numerical 
EN coratron of (2-4) must be considered. 

The errors incurred in a single integration step are 
of two types: 

l. The local truncation error or discretization 
error - the error introduced by the approximation of the 
meererential equation by a difference equation. 

2. Errors due to the deviation of the numerical solu- 
meen from the exact theoretical solution of the difference 
equation. Included in this class are round-off errors, due 
to the inability of evaluating real numbers with infinite 
precision with the use of computer En computers usually 
operate оп fixed word length), and errors which are incurred 
if the difference equation is implicit and is not solved 
exactly at each step. 

If a multi-step method is used, an additional source of 
error results from the use от ап m method (usually 


a single-step method, e.g., Runge Kutta method iPr omen A 


the needed starting values for Тһе тесер ПЕ here 
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A. PROPAGATION OF ERROR IN A ONE SITERE TOD 
The propagation of error in а опе-ссер пегпей сат ре 
analyzed by studying the particular representative one stop 


method of Euler (2-7). 
ere Е р) (2-36) 


This can be done by determining the relation of the error 
Fasten nz] to the error at step n. To do this let ne 
denote the true solution of the initial value problem 

y'(x) = f(x,y) with y (xj) E Then the total accumulated 


Solution error e, at step nis defined by 


(QS 


The numerical values computed by Euler's algorithm (2-36) 


satisfy the relation 


NEED M ба пре Е (2-58) 


where а denotes the round-off errors resulting from 
Ewaiuating (2-36). Similarly the true solution values 
satisfy the relation 

Y 


КО (2-59) 


= ae + DC АКО) + Пере 


where ШЕ denotes the local truncation error in (2-36). 


Subtracting (2-38) from (2-39) yields 


Үні Yne ті ПЕС Se eee 


(2-40) 
where Е +] = TEST * Roel’ Inserting (2-37) in (2-40) yields 
L4] - Y +h [Е (хр, ) - ООС] + ЕТТЕ (2-41) 
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Applying the mean value theorem to the term- im s p 
bracket of (2-41), this relation between successive errors 


can be written as 
E E a ааа (2-42) 


where t, denotes 3f/3y and a lies between Да апа D If 

|£ (x,y) | « C and |E,,,| < E where C and E are both positive 
constants, the error expression above can be replaced by a 
related first order difference equation (Henrici Ref. 5 
wives an excellent treatment of the solution of difference 


equations). 


С ee EM: (2-43) 


Under the stated assumptions it follows from (2-42) that 


ЛІС ІТ ЕТ ЕЕ 


Now if |E,| « e, it follows that 


| E See NI C E 


TEST 
Шбететоте, from (2-43) 

ОРЛИ В е (2-44) 
1t follows then from (2-44) that 


561 < eo > [241 < e, > ....where the symbol > means 


У implies that. 


That is, the propagated error is bounded by the solution of 
the related first-order difference equation, Now since 


МЕРУ. 750 = 08 the condition спас 8 еза рои ленти 


fled by setting РЕТ 0, which is the initial сола олш 


the difference equation. 
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Defining G = 1+hC, the difference едпастол 2 IE a 


be written in the form 


e = бе * E with initial condition e ^ 0 
ПЕЛ n O 


The solution of this tirst-order ditierence Equation еи 


be found by successive substitution as follows 


е, = Ge, + Ë = EÈ 
е, = бе, EP a (GE 
= Е 2 
е: = бе. + Е = (6 +0+1)Е 


е = бе , * E = (601+ rs ,,,+ G+1)E 


The solution e, can be seen to be a geometric progression 


Ane. as such can be written as 





E (2-45) 


го Томо then that the propagated error in Euler 5 alpor- 
ithm (2-4) is bounded by the expression 


rs 1 1 


Pea as a E 


We shall illustrate this estimate by a numerical example. 


Given the initial value problem 
y EOS Oz 


we apply Euler's formula (2-36). It 1s required tomen) mes 
the error at XJ = 1 as a function of h assuming that the 
Не ОБЕ Пе Рог E E 1077 for single precisiongC onn 


precision yields R < 107) in the IBM 360/67 computer. 


S 
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Computing 


С = max Е, (х,у) | 


yields 
С = 1 
The truncation error associated with (2-36) is given by 


Le 
T - 5 ћу (5) пете оиа 


EX 
n+1 n = 


n n+1 


Knowledge of the second derivative y"(x) is needed to eval- 


Mate this. Thus 


1! = d t Зи а 
Е OS a a MXN 
Y Z : 
Since i. = 0 апа > = -], this reduces to 
VO) Te ye 


Using (2-45), one obtains 


A eel 


2 
n= 5; h^ max |y| 


0<х<1 


To evaluate max|y|, it is seen that the analytical solution 


of the initial value problem is 


Then for x = 1] 
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It is clear that max |y| is attained at x =O m s 
nas I 


Substituting into the original expression yic ids 


2 n 
е, 5 [1077 + [em (2-46) 


Rewritine this equation yields 


-7 
e, Š = + 2 Gea eae (2-47) 





Thus the propagated error incurred in stepping forward from 
x = 0 to x = 1, as a function of h, is bounded by the above 
expression. 

The revised form of the propagated error bound shows the 
mir duence the various errors have in the numerical solution 
Ut the ODE. 


Rewriting (2-47) in the form 





B Raed ; n 
y > + bi ((1+h) -1) (2-48) 


where R +1 is still the round-off error 


D is the truncation error reduced ру 0 (t 


and ignoring the contents of the second parenthesis, the 
equation suggests that as h » 0, the truncation error 
decreases toward zero whereas the round-off error increases 
toward infinity. То minimize the error, it is clear that h 
must be chosen so that the truncation error and the round- 


Off error have equal orders of magnitude Алиш аш 
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experiment was conducted to verify this таба РПК ЕЕ 
can be seen that the accumulated error is not simply equal 
to the sum of the local truncation error and round Исая 
but must be computed as given by (2-48) which is the solu- 
tron of the first order difference equation corresponding to 


9-36). 


B. PROPAGATION OF ERROR IN THE MULTI-STEP METHOD 

As shown in equation (2-21) the most general representa 
tion of the Predictor-corrector methods are of the form 
G. y де АН 


1°п-1 3 
0 Т 


Ho 
пи У = 


К 
LO Bi (2-49) 


Yn+1 ` -1 тел 


However by confining the study to a form represented by 


ит = (2-50) 


Simplicity is obtained without loss of generality. It can 
be verified that most of the formulas in Section III are in 
mie form (2-50). For example, the саѕе т = 0 ала 8 1 = 0 


the Adams-Bashforth formulas are obtained as 
У +] - Yn + h 


while in the case r = 0 and 8 1 #0 the Adams-Moulton forms 
result. 

The propagation of error in the malus aman method can 
then be analyzed by studying equation (2-50). Following the 


same definition as in the previous section with regard to 
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Шаты Ls ss на 


, R 1 Е Се апа Е. 


ntl e n+l? ieee lee 


Том Пара 


n 


repetition can be avoided. 
The solution values generated by (2-50) satisfy the 


me lation 


S) 


K 
Уп+] ^ Yn-r * ВО +1>Yp+1? zo о Ват ОК 


1 


- В ЕСІП 


n+1” 


Similarly the solution values Des satisiy the relation 


K 
HO Ë or Ы hB (х,у, үү) + "a BO -4>Yp-3? 
+ Т +1. (2-52) 
Subtracting (2-51) from (2-52) yields the equation 
ee EIL NM. и ҺВ 4 (£5, 1, Yn, 1) 
K 
RE En; 
- f(x, 1) * Dd + R +1 (2-53) 


By application of the mean value theorem to (2-53) and using 


mae definition of 2 and E one Obtains 


+1 
ОЕЕО Е Оер (2-54) 
К = 
m > Ва ¡IVY OY 4) * Винт 


0 
Using the definition of C, E, and Y (2-54) reduces to 


1 


K 
Ener = Inn Beth EB; ¿Ot E 
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OT 


Ë ает C= + dpi 


-1°п+1 n-r ЕЕЕ (2-55) 


n+1 


he related difference equation to (2 55) сапе ЕСИ 


К 
ел-10-Һ18 110) е бг E (2-56) 


ШЕ |7 < е (i = 0,1,...K) where К > r it follows that 


"cn 


n-i 
| К 

| а О ne e ЕНЕР ЕЕ 

Thus 
K 

c S с le. n E 
Erom (2-56) it follows that 

IO ТІР 
If hC[B8. ,] « 1, then 

p (2-57) 


That is, if the magnitude of the propagated error is dom- 
inated by the solution of the related difference equation for 


K*1 successive steps, namely |Z;| < e; (i = 0,1,...K) then 


1 
by induction, the same is true for all successive integer 
 Іпес ой 1. Therefore, a bound for the propagated error in 
the multi-step method can be determined by obtaining a 


Solution of the related linear, inhomogeneous Iene eE 


Equation of (2.55). 
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As with ODE, (2-56) has homogeneous and particular 


sOlutions such that 


e = e tae 
n nH np 


The homogeneous solution is given by 


IU COR K 


е, = (1-hC|B i]|)u U 


E he | 


A еу |! b (2-58) 


The particular solution can be obtained by assuming that 


Em - A, then substituting іп (2-56) 


(1-h|8.1[C) C232 - C2) -hC CO L| 85| *|84] е. 


K 
A[hC[Bg ,| * hC x [8.]] = E 
E 
K 
МС 2 |8.| = E 
т il 
1=-] 
К Е 
Ese ы 
WC 7 В. 
1=-1 Í 


The general solution of the difference equation of (2-56) 
can be written as 
Ж п п п Е _ 
е а uj + duzt... + diu a o (2-59) 
he ENG 


where 41,9,,... де] are determined by the initial starting 


values and where м; (1=1,2,...К+1) are the roots of the 


characteristic equation 


A а [uke |e, fuk ја. + ву = 0 
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In what follows the assumption is made that the roots are 
real and distinct (complex roots and multiple roots are 
discussed in detail in Ref. 4). 

It may be observed that the homogeneous part ot the 
miaracteristic equation for the accumulate dic rro) ти 
exactly the same as in the characteristic equation for the 
original difference equation of (2-50). The significance 
of this point will be used in the analysis of numerical 
stability following a numerical example. 

To illustrate the previous analysis consider the case 


where r = 0 and K = 0 іп (2-50) which yields 


Ерера (2-60) 


го" 


ШІ15 15 the modified Euler formula wherein Це is given by 
the initial condition and the needed additional point has 
meen supplied by a predictor formula. 

Assuming that E, Qo y) - -C and that Е 81 = E, the related 


difference equation yields 
h TS 
(1 2 ЕЗ Uy - 1 - h(-C) DE 0 
which has root 


hC 
р (2-61) 
(1 + >) 


Weing (2-59) the difference equation solution zen 
ne a" 
CONDE 
SPESE ао |e — 1 
(1 + — PA 20 
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Analyzing the вест шыр di as constant, 1t remains to 


show that if 


hC 


е1 
2 


еп the total error е, decreases with n and convergence is 


assured. Furthermore if this condition is met then 


16 
_ 2 1 (26) 
— m Кс - 


ШІ (Пе solution of the difference equation decreases with 
increasing n. But the expression on the left of (2-62) 165 
exactly the root of the characteristic equation of (2-60); 
thus it 1s clear that if uj < 1, then the solution of the 
œ terence equation decreases with n. This idea can be 
Ewtended to the general characteristic equation (2-50) whose 


Solution may be written in the form 


n 


к+1ЧК+1 (2-65) 


DS Eo + dE к 
ШЕ Паз Бееп shown earlier that the total solution error 
Ez e +) is bounded by the solution of the related 
difference equation. In turn it was observed that en depends 
ea the solution of the characteristic equation of the dif- 
ference equation. Therefore the problem of numerical sta- 
bility reduces to the analysis of the roots of the 


characteristic equation of the related difference equation 


Or (2-63). 


Si 





C. STABILITY ANALYSIS ОЕ Я ОЕ ЕР 


Given an ODE 
y! = Ay with ye ya (2-64) 


Assuming that A is constant, the analytic solution is given 
by 
A(x -x ) 
SEE A 
Using Euler formula (2-7), the related difference equation 
! А 
with Yn = Ayp 15 


X SU (1+Ah)y, = 0 
whose characteristic root is 
uj = (1+Ah) 


Therefore the solution is 


ща = (1+АҺ) 


Since there is only one root, (2-63) reduces to 


_ n _ n 
y, = du] = d, (1+Ah) 


where di EUM Therefore, 


HA 
n О 


= yg (1*Ah)^ = у„(1+АһҺ) А (2-65) 





Yn 


We recall from calculus that 


X x/h 


e = Кү 
270 
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and using the fact that A is constant (so ап КАЛ тиан 


no 0), 1:5 сап Пе есен that 


А(х -х ) A(X -х ) 
A AS 
в -0 an 
It follows that (2-65) reduces to 


Аа) 


EK Uc 


The method is then (roughly speaking) stable, which means 
that a sequence of approximations a will have the true 
solution as its limit if the corresponding values of h have 
zero as their limit. This analysis is quite simple since 

no extraneous roots represented by По due are intro- 
duced. The case where these roots are introduced will be 
w cussed thoroughly in the next section. The only error 
introduced is represented by the particular solution of 


(2-59), whose behavior has been thoroughly studied in 
meetion V-A as a function of h. 
ВЕ STABILITY ANALYSIS OF A MULTI-STEP METHOD 
Given the same problem 
г 1 = 

у Ay with y(x,) = Y, 
mme corrector equation of Milne's P-C set (P-C-II) is used 
Ez. 

Yn«1 7 Yn-1 * 3 Yn-1* Yn * Yn+1 СЕ 


Since 
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Y c 


n n 


si и PY +1 


(7-06) has the related difference equation in сес ешр 


У a 48]. y ЕЧ ЛО 7 [ + > - 0. 


The characteristic equation yields 


2 Ah 4Ah АЋ | _ | 
uf [а - AB] му [BE] - [a +B] =o. (2-67) 


with roots 


form 


_ 2Ah #/ 9 + 3An? 


SE АВ 


Assuming that h is sufficiently small, the binomial theorem 
can be used to expand the square root and then simplifying, 


one obtains 


uj = J q Ah + ОВ 
u, s -1 * E Ah * 0(h) 
2 3 
Since n = ғал nk the solution of (2-67) can be written 


as ) 
Бе 
(XD T 2 


y = d} (1+Ah+0(h)) ™ ° + а„(-1)7(1- - Ah + 0(h)) 


As h tends to zero, this solution tends to 


) 


А = 
(ХП Xo) О (2-68) 


у = de + sige)” е 


1 
FA(x, x 
1 ; 
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Since di 1s given by the initial condition sand d, 1S assumed 
to have been supplied by an auxiliary one step method, (2-68) 
can be completely determined. It is clear that there is an 
extraneous root introduced as a result of the use of a E 
order difference equation to represent a first-order dif- 
ferential equation. The root U, 1s then called spurious, 
parasitic, or extraneous and has no relation to the exact 
solution of the differential equation but, nevertheless, is 
unavoidable. This is clearly seen if it is assumed that 
d, - 0, in which case (2-68) reduces to 
АХ 

p o? 
mas is indeed the true solution of the ODE. But in general 
4, # 0, so the behavior of this extraneous root remains to 
be studied, since it affects the total solution. This could 


be done proving that, if A > 0, the parasitic solution 


n n "(X Xo) 
u, = (-1) "ë 
meds to zero exponentially as X increases. Then the 


effect of the spurious solution becomes negligible as Х 


increases while the principal solution 


m IU 
1 
tends toward the true solution values У ОЈ However 1f 


A <0 the spurious solution increases exponentially in magni- 
tude and alternates in sign as the step-by-step calculation 


peegresses, while the principal root и tends torrone: 


SS 





Moreover, at each stage of the calculation попа ки ен не 
will introduce new spurious terms of the same Буре Пе 
these extraneous roots have no relation to the exact solu- 
tion, as they become large, they will dominate the results, 
thus invalidating the numerical solution. This situation is 
meterred to as numerical instability. 

The same analysis can be generalized by going back to 


equation (2-63) 


n 


" n n E 
c diui + du, Бра а l: атча] (2-69) 
is clear then that the principal root is du and 
m n n 
4,1, dzuz,. Фу Uns] ате extraneous roots that are intro- 


@uced as a result of representing a first order differential 
equation by a cea 2 order difference equation. By analyz- 


ing the principal root u, which has the form 


1 


u ТЕЕ а 3.50 (T) 


i 
ет 


cano oh) 


u 
1 
more precise requirements about numerical stability can be 


determined. Observe that un approximates gn as was already 


shown. 
Comparing (2-69) to the exact соттоп оу = ҮЕ 
nAh 
дасе с › 
it is clear that as long as [u,| > [|u;], i 9 2,5,...K*1, 


n 
when n'increases, the extraneous solution represented by Us, 
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will become small when compared to the principal solution 
йуу. In this situation the numerical solution is expected 
to be stable. However, if |u,| > |u,| for any i = 2,3,...K+1, 
u; will become large with respect to u] and the numerical 
solution component corresponding to й w111 dominate the 
solution. This situation is referred to as numerical 
instability. 

Intuitively then, the numerical solution will be valid 
ШЕ Ци. | > Е 1 = 2,3,...К+1. However, recall from the 
рисутоцз analysis of (2-59) that the characteristic equation 


8-59) for the accumulated error e 1s exactly the same 


ntl? 
ие characteristic equation (2-63) for the original 
we rerence equation of (2-50). For a valid numerical solu- 


tion, it has been shown that the total error е, must not 


meow With n. From (2-59), it was observed that 


n n 
e = а. џи. + Ди. +... + d 


п 
nH T 22 he eer 


and is equivalent to the condition that Ju; | < 1, 1525955 ни 
It is then possible to define the stability of a method in 


the following way. 


(1) ABSOLUTELY STABLE if |u,| < 1, i - 1,2,...K*1 


Ші) RELATIVELY STABLE i mu T um s mu 


— 


Absolute stability does not imply relative stability. In 
other words, a numerical solution may have Ju, | sole 
i= 1,2,...K+l, but [м | < |. | < 1, i = 2,3,...Kt+l. 


Summarizing then, it can be seen that in Che TODE 
Vy 


2 





absolute stability and relative stability sare Бол os 
considerations if A < 0 since the exact solution 15 decreas- 
ing with x. If A> 0, the exact solution is growing with 

X,» SO that relative stability is the important consideration. 
In other words, the numerical solution is valid as long as no 
component of из dominates the component corresponding to the 


principal root. 
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VI. CONVERGENCE OP THE CORRECTOR IN ITHE P CENE THOE 


The term convergence has been mentioned earlięr in the 
study of the propagation of the error as a function of h. An 
Malysis of the convergence properties of correctors reveals 
Anat this convergence indeed depends on the step size h. A 
Proper selection of h must therefore be made to ensure the 


convergence of the corrector. 


h 


Let у diee cie WE jt approximation of the solution 


1 


value at x Then the corrector formula can be written 


ntl’ 
in the form with 8 1 Z O: 


Js Lum 
ЕН B % 


I M > 


BER RE BE 


-i’Yn-i) 
(2-70) 


J c 


1= 0 


where 1 is the iteration index. Note that y - УР у апа 


Ax | | | | | 
yi SAS CO re Cto OE предани тие у) which in türn is 


+1 


ШІ corrector оЕ (1-1) iteration. If the sequence of suc- 


А š я + 
cessive approximations у) 


O generated by iterating the 


* * 


corrector, converges to a limit, denoted by y,,,, then У +] 


satisfies the relation 


* 


K 
* 
pd Е ner i ВЕЕ 5 


Ва Ес И И 
(2-71? 


i=0 
Eubtracting (2-71) from (2-70) yields 
Je OS 


q В ) n $ 
Yarl Ons ESA eee Е (Кт? дет? 


which by the mean value theorem can be written 
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Ј+1 2 M СЕ re 
Yn«] Ужа А8 Е one 


j+1 


= * 
where Messi lies between e and Ут 


PE ES < С, where 


к 
Ms a positive constant, inztheznerehberbocczor (c EDDIE 


then 


41 Я : A 
Пе SIDE ES E 


It follows by induction that 


D 


: * : 
j*1 
Dc Ё У 5 ІҺ|В 1161 


О * 
Pru | Xn 
Therefore, the iteration of the corrector will converge to 
the limit E provided that 

ШІСІ (2-72) 


With this expression the conditions for convergence of the 


various corrector formulas considered can be written down 


immediately: 
8 Condition for 
Correc or = Convergence 
Euler | 2 hC < 2 
Milne 1/3 Бе 515 
Hamming 3/8 HC 87/3 
Third-order Adams 5/12 hC 251275 
Fourth-order Adams 2/8 hc < 8/3 


Both the Nystrom and second order Adams P-C sets used the 
Euler corrector while Hermite used the Milne corrector and 


thus are included in the above expressions. 
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An important point that must be brought out here is that 
convergence will only assure that the sequence of iterations 
pea 


Ya+1 Will converge to some definite value y but not 


П 
mecessarily to the true solution y (x,)- Stability is _ still 
the yardstick of a valid numerical solution though (2-72) 
provides a good rule of thumb to follow in the proper selec- 
tion of h. Indeed, later it will be shown experimentally 
that absolute stability can not be attained beyond the bound 
for h given by (2-72) Since convergence is a necessary but 


not a sufficient condition for stability o£ a numerical 


lution (cf. Ralston [Ref. 6]). 
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VII. INFLUENCE OF IHE РЕОРВАСАТЕРИЕРОР 


In order to test the idea presented in the previous sec- 
tion dealing with the propagation of error, the numerical 


example given in Section V-A was run on an IBM 360/67 compu- 


ter where the round-off err bound is по using single 
precision (SP) and J їп double precision (DP). Кеса11 
that the ODE used was 

y' = -y with y(0) = 1 


and the final error expression using the Euler predictor 
formula at x = 1 15 

h 
2 


e = |—— + 


(1th)” - 1 022-733) 





ТО ОАЕ ТЕЕ още 


A series of h values from h = 10 
numerical values of e, Were computed both in single precision 
and double precision. Table 1 shows the actual numerical 
MENS” obtained. These data show that for large h to h = J 
the SP and DP results are essentially identical, but as h 
decreases the round-off errors tend to increase while the 
truncation error tends to zero. The effect of the round-off 
error in DP remains negligible since, in double precision, 
14 digit accuracy is obtained. These results confirm the 
idea that the errors decrease with decreasing values of h 


ЕО a certain point, beyond which the errors increase. A 


plot of the data obtained is shown in Figure 2. 
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TABLE 1 


INFLUENCE OF PROPAGATED ERROR (2) 
FOR y'=-y AT x=1.0 


х 
LÉSp 
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h 
Figure 2. Influence of Propagavead@nr rom 
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VIII. STABILITY BOUNDS O TE 5 


Prior discussion about the stability of the predictor- 
corrector methods relates only to the limiting properties 
of the roots of the characteristic equation as the interval 
of integration approaches zero. In many applications, 
however, additional information about the actual size of h 
1s needed to ensure stability of a method. The number of 
results discussed in published papers on this subject is 
tensive. But the interesting point observed is that 
although they all started with the analysis of the roots of 
the characteristic equations, they varied quite extensively 
Ши пе mode of the predictor-corrector applications amd in 
the number of P-C sets considered. То be more specific, 
some published papers are cited. The real negative stability 
expression used was hÀ where À = t, in a single ODE. Chase 
Metz. 7] analyzed the Milne P-C mode without iteration. The 
The real negative stability bounds found were -0.8« hA «-0.5, 
while on the other hand using the same P-C mode, but now 
iterated to convergence, resulted in numerical instabilities 
for all negative А. Hamming [Ref. 4] used three modes of 
his P-C set. First, iterated to convergence, the resulting 
stability bounds were -0.5 « hA « 0; second, with truncation 
error modification but without iteration the results were 
EN 95 < hÀ < 0; third, still with truncation error molts 


tion but with only two iterations, the results were 
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-0.9 « hà « 0. Lapidus and Seinfeld [Ке TINA 
апа РЕ (СЕ)? methods, where s is the iteration index, pub: 
lished the real negative stability bounds for all available 
predictor-corrector methods. ретт асо сос ае 
different from the algorithm used in this paper in that a 
final derivative evaluatlon is computed before terminating 
the computation with or without iteration. The symbol PECE 
от РЕ (СЕ)? is used to denote their algorithms. The algorithm 
Bensidered in this paper does not require a final derivatıve 
evaluation before termination with or without iteration, and 
hence can be denoted by PEG 

To establish the real negative stability bounds for 
ВР С sets considered herein, the experimental procedure 
used by Chase to obtain the real negative stability bounds 


for hA will be followed. Chase used the test ODE 
у' = 100 - 100у with y(0) = 0 


Using different values of h, he analyzed the behavior of the 
ШОТ until actual instability occurred. The P-C mode will 
edif ferent in the sense that a standard application involv- 
ing two iterations, PEO. will be used for all numerical 
methods considered. The choice of just two iterations was 
incluenced LIinste Dy the publashedepaper of Purik ndier enen 
[Ref. 8]. They conducted an experiment using the PEGI 

mode applied to the Adams' methods only. After analyzing 
several numerical results, they concluded that the best 
method is S = 2, based on the cost of computation, average 


accuracy, and stability. Second, the choice ОТ Устае 


(6 





was influenced by the analysis of the corrector in terms of 
minimum truncation error. Consider Euler's P Cker UM 


which has the form 


Predictor: AA A ћу, 


Corrector: ук: 2 [ "ii 1 


/n*1 n n 


as applied to the ODE y' = Ху. 
Ignoring the corrector, the single characteristic root 


is obtained from the predictor as follows: 
Jn*l ^ 7n ^ A pi. 


Ul 7 (1 + hÀ) = O 


uj = Т thi. 


Now apply the corrector once. Substituting the predictor 


into the corrector yields 
= (а + ћу + 0907 
Yael 5 ( us Yn 


The single characteristic root is 


= 1 + hÀ + 


s ШЕ 
1 2 


Now apply the corrector twice obtaining 











ZZ 
u h мод 
еи SS 
ни ред A~h 
Ур = CL t ha + + 1 І Уа 


mine characteristic root is 
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и] = Ji AA 





2,2 h 


+ 


А 
4 





Continuing this process yields these results 

















O COrrector: u, = l + hì 
20 
l corrector: чу = 1 + ҺА + — 
2 5 
Сове Ове и В н ҺА ҺА 
1 2 4 
2-2 Se 4 
3 correctors: и] = Да a + mA + a 
However, the true solution of y' = ry 1S given by 
22 5. 5 4,4 
¿MA PE п 52 А ~. A 0(һ?) 


and thus the truncation error 


ШЕЛІ eme ру 


т. „РА 
1 
15, first with 0 corrector 
RON 3 
T(x,h) = (1+hA) - Gth o Ch 
Ио 
T(x,h) = DA 0(Һ5) 


By continuing this procedure with the 1,2,3 correctors and 


ignoring higher order terms, there results 


Corrector: IX; h) М 
corrector: = Ween = ee 
2e@orrectors: Т(х,Һ) = 199/12 
3 correctors: T(x,h) = А?Һ?/12 
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It can be seen that the use of the correctorzmor Zee 
seems to yield little advantage in this case. Thus the mode 
Of P-C application to be used will be represented by the 
symbol БЕСТЕ In this mode of application the stability of 
mmc) P-C set depends on both the predictor and corrector 
applications, though more so on the corrector equation. This 
assertion will be clearly illustrated in the case of the 
Milne, Hamming and Nystrom P-C sets. Thus, a priori know- 
ledge of the stability of a method can be obtained by 
analyzing the roots of the characteristic equation of the 


eorrector formula. 


A. EXPERIMENTAL PROCEDURE 

Two general algorithms were developed to implement the 
Ша шистог-соттестот methods in Section IV-E. The first 
Iori thm covers seven of the P-C sets, while the second 


mecgorithm covers the E 


Р-С set, i.e., the Hamming modified 
predictor-corrector set. The flow chart for ALGORITHM 1 15 
расте 164 and that for ALGORITHM 2 is on page 167. How- 
щеш, the general steps for each algorithm will be outlined 
пете. 

Given the necessary starten? values, һе Step size, the 
Gonvergence term and the range of integration, the iterative 


Predictor-corrector method for computing the numerical solu- 


tion а 1s set up as follows: 


ALGORITHM 1: 


Step 1: Compute the predicted vale’ by Ириска 
formula. 
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Step Compute the derivative of the predicted valne: 

Step Compute the corrected value by Che correas 
formula. 

Step Test for convergence. If convergence is attained 
go to step 7, otherwise proceed to thorn pot opa 

Step Compute the derivative of the corrected value. 

Step Compute the new corrected value by the corrector 
formula, using the value computed in repr sas 
the new predicted value. 

осер The result is taken as the desired solution 


ALGORITHM 2: 


point. Advance the integration point by the step 
size. Return to step l to compute ие Еос 
solution point. 


tep 1: Compute the predicted value by the predictor 
formula. 
Step Modify the predicted value by adding the trun- 


cation error* of the predictor formula. 


Step Compute the derivative of the modified predicted 
Value. 

Step Compute the corrected value by the corrector 
formula. 

Step Modify the corrected value by adding the trun- 
cation error“ of the corrector formula: 

Step Test for convergence. Тї сопует у енее не а лыны 
go to step 10; otherwise proceed to the next 
step. 

Step Compute the derivative of the modified corrected 
Matices 

Step Compute the mew corrected valero NCC OM ME 


* 


formula using the valiievas vcomputed necs 
as the new modified predicted value. 


Truncation error as used in the formula is not the 


Original truncation error but the modified truncation term 
expressed as a function of the difference between the cor- 
rected and predicted values. 
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Step. 9: Modify си Ес есес ае 
Step 10: The result obtained is accepted as (lop ок 
solution point. Advance the integrationm Donni 


by the step size. Return to step 1 to compute 
the next solution point. 


In both algorithms the mode of application 1s PEO as can 
EE cadily verified. Іп пас señeme,. 1 the Сеттес ко е 
need only once, then the number of function evaluations needed 
is also one, since it is assumed that the function value for 
the corrector has already been computed; if the corrector is 
used twice then two function evaluations are needed. Thus 
the number of function evaluations is determined by the 
number of iterations. 

The P-C algorithms were coded using the FORTRAN language. 
EN ortran programs for Р-С-Т to P-C-VITI are listed on 
pages 170 to 200. In each program the meanings of symbols 
and parameters are explained through narrative comments. 

To determine whether the stability limits have been 
reached, certain criteria must be followed. Often, when the 
stability bound is approached through increasing step size h, 
the method begins to lose accuracy. An inaccurate, though 
stable, solution may often result in oscillations which 
appear at first to be due to actual numerical instability. 
Actual instability, associated with too large a value of h, 
usually results in an increasing oscillation in the error, 
or simply, the growth of the error in one direction. The 
Particular error behavior depends on the sign of пора 


sıtic root causing the instability; a negative roo ticas 
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oscillations and a positive root causes a miro rn y 


growth. 


B PREDICTED STABILITY CHARACTERISTIC Ор THEE CR PI 

To gain a priori knowledge of the behavior of the sta- 
bility characteristics of the P-C sets, the roots of the 
characteristic equations of the corrector formulas will be 
Studied. It must be noted that the behavior ot the charac- 
teristic roots, though dominating the P-C method used as a 
mer, will be influenced by the behavior of the predictor 
formula. Since the chosen mode of application, DE Í 15 
quite different from ali the other modes of applications 
published in different papers, the stability limits obtained 
Y Several authors could not be used as the stability limits 
of the P-C set considered here. Thus the actual real nega- 
tive stability limits will be determined through numerical 
experiment. However, as noted previously, the analysis of 
the characteristic roots remains the same, and as such the 
published results are used and references are cited. This 
15 done solely to conserve space and to avoid the tedious, 
mepetitious process needed to solve the characteristic 
equations. The procedures have been amply presented and 
mielustrated in detail in the section on numerical stability 
of predictor-corrector methods. 

a) Analysis of the characteristic roots of the со пишеше ол 
formulas. 

The P-C sets with common correctors ма шестој са 


together except for the case of the second order Adams' 
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method which has the same corrector as Euler, and which will 
be included in Adams' type. 
Euler (P-C-1) and Nystrom (P-G IIT) 


The corrector formula is glven by 


ED 


+ 


Dee 
Mie” related difference equation 1s in the form 


TUNE | hC _ 2 
(1 >) т (1 f >) Yn 0, where C Е. (Х,У) 


Me single characteristic гоої 15 


ТШ 
2 

m 
1 hC 
iz 


It can be seen that if C < 0, (1.e., the derivative of/dy is 
Шерастуе), then for absolute stability to occur, |у | < 1 
Mich is satisfied 1f hC/2 < 1. Thus the numerical solution 
will decrease as does the exact solution ТЕ С > 0, ав long 
as hC/2 < 1 then the numerical solution increases as does the 
ВЕС solution. Figure 3 exhibits the behavior of the root 


versus hC as shown by Lapidus and Seinfeld [Ref. 2]. 


Hiline (P-C-II) and Hermite (P G IV): 


The corrector is of the form 


e h : 
Yet Улу з Ор ELEME 


The related difference equation is 
io ДАС hC _ 
| MO Уд ` ЕЗ Yn 7 m E 
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Figure 3. Characteristic Root for Euler Corrector Formula. 
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The roots of the characteristic equation, as derived in 


section V-D, are of the form 


51-801 ШЫ 
1 
3 


= -] + 


u, hC + 0(h) 


which can be written as 


" „hC 
i 

_ _-hC/3 
une 
ШЕ С > 0, um 


| n 
се, ЕЯ < 1. When, however, C < 0, u, decreases as 


behaves like the exact solution and u> dies out 


does the exact solution, but u, increases. Thus the corrector 
is relatively stable but for C > 0, has no real negative 
stability bounds. Figure 4 shows the behavior of the roots 
as given by Chase [Ref. 7]. 

Hamming (P-C-V): 


The corrector is given by 


Sl _ 3h Р _ 
Е Те Рот t; i 2У1 na 


The related difference equation is of the form 


oo| — 
“2 

|| 

о 


Yn+1 Í t 4 8 |^n-1 n-2 


3hC 
еі 





IA Ек 
п 





The characteristic equation is 


8 8 4 


15 k | | m Ë ‚ shc) h 


3nC би 
el: 


Chase [Ref. 7] analyzed the root behavior of this character- 


istic equation. Figure 5 shows the behavior of these 
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Flgure 4. Characteristic Roots for Milme Corrector formulae 
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x - Positive Roots 


N - Negative Roots 


O - Complex Roots 
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Figure 5. Characteristic Roots for Hamming Corrector Розы вай 
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characteristic roots versus hC. Analyzing пе рү арш 
Figure 5, it is evident that the corrector но не ва о Pp r 
22.4 «hC « 0. This corrector has also Zr de zer ner wer: 
relative stability as Chase has shown. In the graph the 
symbol used for positive roots is an x, that for negative 


meets 15 an N, and for the complex roots is an O. 


ENMIAMS-MOULTON CORRECTORS (Second Order (P-C-VI, Third 
ШЕСІ (P C VII), Fourth Order (P-C-VIIT)): 


Since these formulas are all of the same type, they will 
De studied together and extended to the general case of the 
Adams-Moulton types. For the second order, the corrector is 


given by 
Уд 1 ^ Yn * 7 ШЫ ШЫ 
The related difference equation is of the form 
1-39 yaa 01 + 36) у = 0 
meen characteristic equation 


[1 - 25] и - (1 + 25] = 0 


in the limit, ав h > 0, the root becomes 


mor the third order, the corrector is 
у = у + h [5 1 ЩЕ 8 1 Е 1 ] 
ntl n E? Yn+1 Yn И 


Me related difference equation is given Dy 
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5 2 ҺС _ 
[1 - 154 hC] y. - Il + he E EA s 


I 


with characteristic equation 


5 2 


hC 
BN. 12 hC] uq I R: 


hC] u + — 


о 
|| 
С 


Mithe limiting case h +> 0, пре roots are 


и(и-1) = 0 


Рот the fourth order, the corrector is 


I h t t Ях 1 
use] ^ Xue 2% PS esl š 19У п il AE 


The related difference equation 15 


3hC 19 ShC 





[1 - ===] У UE 2c i = 
8 n+1 24 n 24 В p PN 
with characteristic equation 
3 3 19 2 Shc ИС а 
[1 я 8 hC]u E [1 + 24 ВС] ци + 24 li D = 0 


M the limiting case as h > 0, the roots are 


u* (u-1) = 0 


Thus, it can be seen that, in each case, the dominant root 
15 1 while all the other roots (the parasitic ones) lie at 
the origin when h = 0. In the general q-step Adams-Moulton 


formula the equivalent characteristic equation is 


ne = 0 
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Thus, in every case, the A-M formulas have one root on the 
unit circle with all others at the origin when h = 0. 

As such they would seem to have excellent stability charac- 
teristics no matter what q-step formula is used. When h 
increases from zero, the parasitic roots move toward the 
unit circle. However, significantly large values of h can 
be reached before instability occurs. Crane and Klopfen- 
stein (Ref. 9] analyzed the behavior of the general Adams- 
Moulton formula and showed that it has stability bounds 
МЕР < С < 0. It is interesting to note that, in general, 
as the order of the Adams corrector increases, accuracy 
increases but stability decreases. This is the main reason 
why higher order formulas are not considered (1.e., fifth, 
Sixth, seventh, eighth). The graph of the characteristic 
roots of the general A-M formulas, as presented by Crane and 
Klopfstein, was not reproduced since it would give no 
additional information. 

In the foregoing analysis of the characteristic roots, 
the predicted stability limits served as a guide in the 
choice of the step size h to be used in the numerical 
experiment for obtaining the actual real negative stability 


bounds for each P-C set considered 


в. NUMERICAL RESULTS FOR REAL NEGATIVE STABILITY LIMITS 

In order to obtain actual real negative stability bounds 
See the P-C sets considered in the P (EG mode of application, 
the ODE 


у в -y with y(0)= 1 
was chosen, where fy (x,y) = C = -1 < 0. 
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Each P-C set was run with increasing value moimoi 
size h, in small increments, to determine more precisely the 
actual stability limits. The criteria for actual numerical 
instability set forth previously, were followed in the 
analysis of data obtained. The prior knowledge of the roots 
of the characteristic equation obtained in the previous 
analysis narrowed the choice of starting values for the step 
size h for each method. Each method was run in single 
precision because of the large values of h. 

1. Euler (P-C-I) 

Since the predicted stability of the corrector is 

the range -2 < hC < 0, the starting value chosen for h 
Ша 1.1, then was increased by 0.1 until instability occurred. 
maple 2 shows selected results for different values of h. 
At h = 1.1, it is absolutely stable since the error continues 
to decrease as x increases. The first sign of oscilla- 
tion was observed at h = 1.3 but it is still stable because 
БИЕ error is decreasing. Oscillation becomes evident at 
m= 1.6, but although the solution is inaccurate, it is 
not yet unstable, because there is no uniformly increasing 
trend of oscillation. Actual numerical instability occurred 
at h = 1.9 as is clearly shown by the uniformly increasing 
oscillation. Thus this method is stable within -1.9« hC« 0. 
As compared to the predicted stability region -2.0 « hC « 0 
obtained in [Ref. 2] where in the corrector is iterated to 
covergence, this was expected because of the use of only 


two iterations. 


81 





TABLE 2 


ABSOLUTE ERROR (2) FOR EULER (P-C-1) 


ODE: y' = -y with y(0) = 1 
Single Precision 


[ana T = 
ee > 2 u 
К 2 


8.023270 


2.6 -1.019880 - 
3.9 ne 
a e 
6.5 “КЕТО О 
7.8 a ESTO 
9.1 CODES 
(STABLE) 
ћ = 1.9 
x ЈЕ 
"шый Т 1.9 -6.696-10 
4411072 3.8 5.308:10 ! 
449-107? Su 1.546 
.894:10 2 7.6 3.020 
-1.442:1071 9.5 -6.363 
CS (UNSTABLE) 


(STABLE BUT INACCURATE) 








22 Milne ACA 

A priori knowledge indicates that this methoden Ms 
real negative stability, so the starting value of h was 
chosen as 0.1. Table 3 shows the selected data obtained. 
Even at the starting value at h = 0.1, the buildup of error 
is consistent though small. Then at h = 0.4 oscillations 
become evident and the error tends to increase as x becomes 
large. At h = 0.7 the same increasing error tendency was 
observed. These observations, based on actual results, 
confirmed the predicted instability of this method for C <0. 
Thus it can be concluded that this method has no real 
negative stability bounds, as is evident from the actual 
Aata presented on Table 3. 

s Nystrom (P-C-IIT) 

The selected starting value is close to the value 

Са іп Euler since sthis method has the same corrector. The 
program for Nystrom was run starting at h = 1.0 in increments 
of 0.1 until instability occurred. Table 4 shows the actual 
@ata Obtained. The result for h = 1.1 to 1.3 show absolute 
Stability. For h = 1.5, the error starts oscillating but not 
at a uniformly increasing rate. The solution corresponding 
to this data is inaccurate but stable. This inaccuracy 
continues up to h = 1.7. For h = 1.8 actual numerical 
instability has occurred as shown by the increasing oscilla- 
tion. Thus this method is stable for -1.8 « hC « 0. This 
range of stability is less than that for the Euler method 


though they have the same corrector equation. This confirmed 





TABLE 5 
ABSOLUTE ERROR (Z) FOR MILNE (P-C-II) 


ODE: y'=-y with y(0)=1 
Single Precision 


aan 
S 
(UNSTABLE) 
m 
4. 


(UNSTABLE) 





X ay 

Pee КООШ е 
4.2 gue eng ^ 
5.6 SR 
7.0 1.186°107° 
8.4 1.918:10 7 
9,8 рола Іше 

UNSTABLE 
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TABLE 4 
ABSOLUTE ERROR (Z) FOR NYSTROM (P-C-III) 


ODE: y' = -y with y(0)=1 
Single Precision 


4.488* (STABLE) 


Zi 
(STABLE) 


7.005:1072 


6.538'10 * 


8.581'10 * 


-2.784'10 ° 


1.128'10 | M 


(STABLE BUT INACCURATE) (UNSTABLE) 





бо 





the idea that the stability of the predictor corr шыш 
depends on both the predictor and comrecton ironia 
4. Неше (F CRIN) 

This method has the same corrector as the Milne 
method, hence the choice of starting value is the same, 
Du 0.1, which is then incremented by 0.1 until instability 
Becurred. Table 5 shows the selected actual data obtained. 
From these results it is evident that this method xs unstable 
since for all values of h the error growth is steadily 
increasing. This is expected since it was predicted that 
bue corrector has a dominant role in the stability of this 
method as a P-C set. Thus this method, which uses a Milne 
porrector, followed the instability behavior of the Milne 
corrector, within two iterations. Therefore, it can be 
seen that the Hermite P-C set has no real negative stability 
bounds. 

5. Hamming (P-C-V) 

The chosen starting value was h = 0.5 with increment 
ШОШО ine predicted stability of the corrector 15 2.4 <hC< 0: 
0.6, the method is 


Table 6 shows, that for h 0.5 to h 


|| 
| 


absolutely stable. For h 0.7 to h = 0.8 the method is 
Stable although the numerical solution becomes inaccurate. 
For h = 0.9 the error grows in one direction thus resulting 
in actual instability. Thus the method is stable within 

me range of -0.9 < hC < 0. | сто Ке релге Ее 
ШЕНІ Су ог the Hamming corrector, given by -2.4 < по < 0, tne 


range of stability is greatly reduced. This is due to the 
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TABLE 5 
ABSOLUTE ERROR (F) FOR HERMITE (P-C-IV) 


ODE: y'=-y with y(0)=1 
Оре Рес сло 


(UNSTABLE) 





x SEE 
1.6 eo 
Ж; GROS 
4.8 SEE 
6.4 O 
8.0 BSO 
9.6 1.496°10 * 


(UNSTABLE) 


ov 





ТАВІЕ 6 
ABSOLUTE ERROR (£) FOR HAMMING (P-C-V) 


ODE: y'=-y with y(0)=1 
Single Precision 


4 4 


-3.741'10 


-2.056:10 


4 4 


ОЕ -2.897'10 


1061102 rosa | 


EE | КОЛГО 


62107? 


78107? 


181-1072 


.459*10 ^ | BOUES 


3 


.096* 10. оО 


Fa | ща ООО. 


112171072 | с рана 


‚874° | -1.534°10 > 


(STABLE BUT 
INACCURATE) 


(UNSTABLE) 
.001' 


(STABLE) 
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effect of the Milne-predictor. But recalling the discussion 
of the published results, this actual real negative stability 
bound is exactly the same as that obtained by Hamming using 
two iterations with truncation error modification. This 
coincidence is not without basis since the P (EC) mode has 
exactly the same truncation error modification as that of 
the Hamming method, making the two experimental procedures 
identical. 

6. Second Order Adams (P-C-VI) 

The starting value chosen was h = 0.1 with an 
increment of 0.1. Table 7 shows the actual data obtained. 
For h = 0.1 up to h = 0.3 the method is stable. But for 
h = 0.4 the method shows actual instability as evidenced by 
the uniformly increasing error growth. Thus this method 
15 stable within the range -0.4 < hC < 0. Compared to the 
predicted stability of the corrector, which was the same as 
Ш: (P-C-I), with range of stability -2.0 < hC < 0, it is 
evident that the stability is greatly reduced. Again this 
memld be traced to the high instability of the predictor 
formula used. This is probably the reason why almost all 
published material on the Adams-Moulton method always starts 
mech the third order and higher predictor. Hull and Creemer 
[Ref. 8] showed that this method has the largest error among 
the Adams-type methods. 

7. Third Order Adams (P-C-VII) 
The choice for h starts at h = 0r and i Ене 


incremented by 0.1. Table 8 shows the selected actual data 
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ТАВБЕ 


ABSOLUTE ERROR (X) FOR 5ЕСОМЬВОВОЕВ ВА ЕН 


ODE: y'=-y with y(0)=1 
Single Precision 


1:315: | ШІП 
(STABLE) (STABLE) 


. 6061: 
7505" 
5706066: 
5.060808: 


3.868-10° | 9. 
(STABLE BUT INACCURATE) (UNSTABLE) 
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TABLE 8 
ABSOLUTE ERROR (X) FOR THIRD ORDER ADAMS (P-C-VII) 


ODE: y'=-y with y(0)=1 
Single Precision 


1.600-10 2 


8.272'10 7 


1.832:1073 


4.614: | ZORRO 
4.608- (STABLE) 


(STABLE) 


-4.736'10. 
11765103 
1.569:10. 
1.869°10 


(UNSTABLE) 


(STABLE BUT INACCURATE) 
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Ebtained. For h 0.8 to h = 0.9 absolute stability есе 
For h = 1.0 to h = 1.2 oscillation occurred and the solution 
is inaccurate though stable. For h = 1.3 the error growth 

is increasing showing actual instability. Thus the method 

Ms table within the range -1.3 < hC < 0, which is in 
conformity with the predicted stability of the general Adams- 
method, -1.3 < hC < 0. This also justifies why most Adams- 
Moulton methods start with the use of the third order predic- 
por. 


ees Fourth Order Adams (P-C-VIII) 


Starting value was h = 0.5 with 0.1 пресвете 


ole 9 shows the actual data obtained. From h 0.5 to 


һ = 0.7 the method is absolutely stable. For h 0.9 oscil- 
lation continues as with h = 0.8. This behavior was observed 
üp to h = 1.0 but the numerical solution is still stable 
though inaccurate. For h = 1.1 actual instability is evident 
from the tendency to increasing oscillation of the error. 
Mus this method is stable for -1.1 < hC < 0. 

Summarizing then, the real negative stability bounds of 
the different P-C sets considered are outlined below with 
C = Е (ку) and h = step size. 


Real Negative 


P-C Set Stability Limit 
P-C-I -1.9 < hC < 0 
Р-С-11 Unstable 
P-C-III -1,8 < hC <'0 
P-C-IV Unstable 
P-C-V -0.9 < hC < O 
P-C-VI 0, AAC ЮО 
P-C-VII Е е 
P-C-VIII -1.1 < hC < 0 
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АВЕ 


ABSOLUTE ERROR (2) FOR FOURTH ORDER ADAMS TEE ШЕ 


ODE: y'=-y with y(0)=1 
Single Precison 


2,125:107 


1.488:1074 


sil oe 


3.95.1077 


=. 


зоо [S 


6 


7.9:10. (STABLE) 


ОЛ Р 


des 


(STABLE) 


1.819-10°° 


2.111°107° 


4.93-10^? 


9.0 а ево тај ~ 


(STABLE BUT INACCURATE) 


И 
(UNSTABLE) 
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The significant points noted from these res eee 


1) Within the limits of their stability, the Hamming 
and Fourth Order Adams methods produced the best accuracy. 
11) It is quite interesting to see that although the 

Euler and Nystrom methods have the widest range of real 
negative stability, their results are not quite as accurae 
as the results of the Hamming and the Fourth Order Adams, 


Een the Third Order Adams methods. 


t= NUMERICAL RESULTS FOR RELATIVE STABILITY 

If C» 0, the solution itself and generally also the 
от are increasing exponentially. In this situation 
relative stability is the important consideration. А Р-С 
method is relatively stable if the rate of change of the 
error with respect to a finite range of integration points 
15 less than the rate of change of the true solution with 
respect to the same finite range of integration. By this 
definition it could be seen that it is extremely difficult 
to establish a fixed bound for relative stability. In 
order to have a working knowledge of the relative stability 


of the P-C sets considered the ODE 
y' "y with y(0) = 1 


where C - 245422 > 0 was used. Each P-€ set was run with 
the same series of h values from h = 0.1 to h = 2.0, with 
an increment of 0.1 The range of integration for each 
step s1ze h was x = 0 to x = 10. Table 10 shows the true 


solution values from x = 1 to x = 10.0. Table 10-A to 
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10. 


TABLE 10 


TRUE SOLUTION VALUES FOR 


ODE: 


y'=y with y(0)=1 
comode meu Su 


Yexact 


20. 


54 


148. 
403. 
1096. 


2960 


8102. 


27200455 


ЈЕ 


Шо 
525,9 


085 


2597 


409 
417 
595 


„858 


710 


520 





Table 10-H show the selected actual data obtained нон 
different P-C sets. 

Table 10-A shows that the Euler method is relatively 
stable for the particular range of integration. Comparing 
pne error growth with the solution growth of Table 10, it 
could be seen that the rate of change of error is less than 
што rate of change of the solution. From h - 0.1 to h = 0.5 
the method is quite accurate. However from h = 1.0 to 
ШЕ 7.0 1 becomes quite inaccurate though still relatively 
stable. Table 10-B shows that the Milne method is not only 
A latively stable but accurate from h = 0.1 to h = 2.0. 

From Table 10-C the Nystrom method exhibits relative stabi- 
су from h = 0.1 to h = 1.0 but becomes unstable at h = 2.0. 
This could be determined by simply comparing the magnitude 

8 (Пе error at x = 10.0 апа h = 2.0 with the magnitude of 
mie true solution at the same point, in which case it is 
obvious that the error is larger than the true solution. 

This is another way of looking for relative instability 

Шише it is clear that if the rate of change of the error 

15 greater than the rate of change of the solution with 
respect to n (the number of solution points) then eventually 
[Пе error will be greater than the true solution. Table 10=D 
Shows that the Hermite method is relatively stable and also 
accurate. From Table 10-E it can be seen that Hamming's 
ШҮспосй 15 also stable and accurate. "Table I0-F shows thart 
the second order Adams’ method is stable from h = 0.1 to 

lm = 1.0, but the accuracy is quite diminished av AA 


At h = 2.0 it is unstable. From Table 10-6, the third 
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TABLE 10-A 
ABSOLUTE ERROR (2) FOR EULER (P-C-1) 


ODE: y'=y with y(0)=1 


Single Precision 


-4.819°10 -1.567:10_ 
-2.570'10_ -7.515'10. 
-1.046 - 2.959 
- 3.808 -10.638 
-13.011 -36.261 


4210721 1195555 И. 


. 360 - 136.450 -383.286 а 


19.774 -427.070 12082400 -1420.042 
-60.496 Lolo 17 = 908200 = 


1823828 400602574 -11546-150 О Ü 
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TABLE 10-В 
ABSOLUTE ERROR (X) FOR MILNES E C DIS 


ODE: y'=y with y(0)=1 


Single Precision 


И 


56472.7027 
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TABLE jo 
ABSOLUTE ERROR (£) FOR NYSTROM (P-C-III) 


ODE: y'=y with y(0)=1 


Single Precision 


-2.756°10 
-2. 23110 -5.244*10. 
-1.016 -2.758 
-3.901 -11.304 
. 096° -13.757 -41.718 
. 996 -46.134 -145.106 -215.571 


= 05.900 119 050 - 4852911 = 


ZI EM = NN 215 82.0757 -2584.042 


00.2383 - 14798559 ESO 115 = 
NS Orel Oo o E 5 -1539/57800 -2442:5330 
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IABLE 10 D 
ABSOLUTE ERROR (2) FOR HERMITE R PECI) 


ODE: y'=y with y(0)=1 
Single Precision 


1 


7.092:10. 


ААИ 


2860.671 








TABLE 10=Е 
ABSOLUTE ERROR (5) FOR HAMMING (P-C-V) 


ODE: y't=y with y(0)=1 
Single Precision 


AA 


202: 4.332:1071 


1 


.004: 8.498:10 


ELE 1.043 

. 039 -5.24:10 1 

. 639 10.532 380.836 
Se : 


- 21 38566 3040.449 
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TABLE 10-F 
ABSOLUTE ERROR (£) FOR SECOND ORDERTADANS tar 


ODE: y'=y with y(0)=1 


Single Precision 


-2.735 "10. 


ISO OO 


«I. JT =2 E 

-4.453 2 149599 

-16.013 -45.545 

-54.444 -161. 08 z 20] MT 
=O: = 3А 254 083 554 = 
= 2147696 - 5005954 -1812. 745 -2524.0472 
-06% 800 IAS SÓ 2:5 909051957 - 
202 378 = ои 92 -18684.180 -—— 1 06025 Ü 








TABLE 10-6 
ABSOLUTE ERROR (í) FOR THIRD ORDER ADAMS (I GQ iS 


ODE: y'=y with y(0)=1 


Single Precision 


2 d sese 


Enn 
‚569. 
.062: 
ІШЕ 


. 524: : 11.147 


1412: 27.090 39,645 


507 SUS 132.954 2172016 
. 089 209.992 ALS = 


82828 5 #549. 750 2409.128 
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TABLE 10 H 
ABSOLUTE ERROR (2)JSFOR FOURTHS ORDENADA MS ТР БРИ 


ODE: y'=y with y(0)=1 
Single Precision 


7 


9.992:*10 


-7.852:10 ` 


=o Od 

220350 
- 177926 2819859 202 2008 
2509578 "25 5220 > 


= 709812 930.164 1748.47 
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order Adams' method showed accuracy and stability from 


nos 0.1 to h 


2.0, Finally, Table lO0-HosShowsmbh o ME 
Fourth order Adams' method is quite stable and accurate 
írom h= 0.l to h= 2.0. 

The most interesting points observed from these actual 
mesults are: 

i) Though the Milne and Hermite methods have no real 
negative stability limits, they are quite accurate for 
С > 0, In fact the Hermite method produced the least start- 
ie error, and the Milne's method provided ehe second least 
Siearting error, and both maintained accuracy up to h = 2.0. 

mie in contrast, the Euler and Nystrom methods both 
have the widest range of real negative stability limits but 
Showed inaccurate results starting at h = 1.0. In fact the 
Nystrom method is unstable at h = 2.0. 

111) The third and fourth order Adams' and Hamming methods 
Showed wide range of relative stability and provided accurate 
results. The second order Adams' method again showed inac- 
curate results at h = 1.0 and relative instability at h = 2.0. 

iv) From h = 0.1 to h = 1.0, the methods in terms of 
accuracy, rank as follows: Milne ranks first, Hermite second, 
then Hamming, Fourth Order Adams, and Third Order Adams in 


that order. 
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Tk... NUMERICAL EXPERIMENTS 


In order to test the performance- of the different P- C 
sets considered, a collection of test ODEs with many unusual 
ema interesting features (i.e., singularities, discontinui- 
ШОС Infinite derivatives, oscillating derivatives, etc.) 
were selected. Several of such ODEs were presented by Hull 
and Creemer [Ref. 8] and Lapidus and Seinfeld [Ref. 2]. 
Table 11 lists the test ODEs considered. To simplify nota- 
mon, the ODE number will be used to specify the differential 
equation. For example ODE I is equivalent to the ODE 
DES -y + 10sin3x with y(0) = -3, whose analytic solution is 
BEES = sin3x - 3cos3x. 

From the previous section on stability bounds of P-C 
sets it was observed that the different methods showed good 
meeuracy up to h = 0.5 on both experimental ODEs, y' = -y 
mec y' = y, considered. It was further observed from pre- 
wus numerical results that all the P-C sets exhibited 
mood stability behavior within the range of h up to h = 0.5 
except for the Milne and Hermite methods, in the case of 
real negative stability. But these shortcomings of the Milne 
and Hermite methods were compensated by the fact that they 
produced good accuracy and wide range of relative stability. 
This analysis is needed in the sense that the test ODEs con- 
mered exhibit both cases of Еу (х,у) < 0 and f (x,y) > Чу 


= 276 = - 4 => =A =1 


By using values of h= a E ON ZN 
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where the largest value of h = 0.5, each method will Бе 
expected to perform quite well. Thus a comparative analys ik 
of their accuracy and computing time woulda be meaningful. 

In dealing with accuracy, the effects of the propagation 

of errors would be strongly felt, as has already been shown 
in the previous analysis of error propagation. While 
computing time takes into its fold the number and complexity 
of function evaluations. 

АШ кше Мена аре Е ЕБЕ уеге mui coches tT ENS 
ie Series of h values as previously mentioned. The Fortran 
programs used were the same as those used for stability 
analysis, but the precision used was Double Precision (14 
digit precision for IBM 360/67 computer) to minimize the 
ЕЕС of rounding error as h decreases. Each test ODE will 
be studied individually, with the objective of providing 
useful and important recommendations on which P-C set 


performs best for the particular class of problems. 


АО ODE I 

Table 12 shows the exact solution from x = 1.0 to 
EN U.0. From Table 12-A the results show that the Euler 
method provides the best accuracy with 7.65 secs in computing 
time with the Nystrom's ranking second in accuracy though 
puth less computing time, 7.46 secs. The Milne and the 
Hermite methods are not considered since they both showed 
instability in their numerical solution values. This could 
easily be seen by looking at Table 12-B which shows the 


behavior of the errors for these methods. The errors for 


JO 7 
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TABLE 12 


TRUE SOLUTION VALUES FOR ODE 1 
Double Precision 


2 Y exact 
1.0 3.111097 
20 А. 
5.0 3.145509 
4.0 - 3.068134 
5.0 О ОЗ 
6.0 32.034957 
uo 2.479843 
8.0 “2175115 
9.0 127092792 

10.0 -1]. 450765 


109 





TABER PZA 
ОВЕ I AT x = 10 0 


ABSOLUTE. ERROR FOR P GC | ° P- CS P G Ill an C BEI 


Double Precision 


“HE 
-9.83522' ооо" 
93415: -3.40492°10°° 
METUS. Ела seo osa 
.92664“ о, 
(522578 Eo sass 
35324: 00786: -6.039341 
TOO .38212' -37.089294 
46 E 





я 
CPTS means computing time in seconds. 


«я 
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TABLE 12:15 


ABSOLUTE ERROR FOR О Р=С=1, P-G TIT SP C lll o Ban ЕЕ 


pouble Precision 
ODEI; h= 0% 





Пе 





the Milne and Hermite methods increased in one direction 
eventually growing to magnitudes greater than the true 
Solution at x = 10.0 (1.e., for Milne's error Бе 
and for Hermite's error -37.089 >1.450), thus these methods 
are not suited for solving ODE І. The Euler and Nystrom 
methods both started at accuracy approximately equal to 
fifth decimal places and ended up at h = 0.5 with accuracy 
up to the first decimal place. But though their accuracy 
ES quite restrained as h increases, nevertheless the behavior 
Se the error tends to decrease as the range of integration 
Mseincreased. From this group then the Euler method is 
ШІСІ 265%: choice for ODE I. 

table 12-С the results for the other four methods ean 
listed. All these methods showed stability for all values 
of h used, but Hamming's produced the best accuracy and 


the best computing time. Hamming's method started at Ыш 


оо ие 


 Ситасу at 1/h = 128 and ended with 10 
Fourth Order Adams' comes next in both accuracy and computing 
mes. then the Third Order Adams', and the Second Order 
Adams' came in that order. 

Comparing these last four methods with the Euler method, 
only Adams' second order method is inferior. Therefore, for 
solving ODEs that belong to the class of ODE I the Hamming 
method is highly recommended in terms of accuracy and least 


post in computer time. The Fourth Order Adams’ 15 the next 


choice in the absence of the Hamming method. 





TABLES? -C 
ODE I at x = 10.0 


ABSOLUTE. ERROR (59) БОҚ Р-С У, Р-с-үгарас SA 


Double Precision 


10 |-9. 86551. -9.60790* 1.86070-10^ 


3.0381:10. 


8 


.12626-10 „95790: -7.81167: 18596: 


= 


.68494:10 E17 005 155509, .04497° 


‚01544°10° 


212959: . 49889" .91417* 
2997024 B5 Зоро що ‚464035 
2591507 051707 .00046° S0055 
#5990607 .586413 SUI OS гащи аа 


27.5 -90 25. -99 








Be ODE II 

The true solution values are given in Table Rom 
x= 1.0 to x = 10.0. From the data presented in Tabi Ek ai 
it is clear that the most important consideration is the 
Шер Size h. Рок -h-= 0.5 the foun methods identically 
yield errors greater than the true solution at the same 
point. In solving problems of this type then h must be 
chosen to be smaller than 0.25 for the numerical solution 
to be valid using the Euler, Milne, and Nystrom methods. 
The Hermite method obviously is not worth considering, 
because of its inaccurate results for all values of h. 
For h < 0.25 the Milne method gives the best numerical 
solution and the least computing time, followed by the Euler 
then the Nystrom methods. 

Table 13-B shows that to have a vlid numerical solution 
for any of the other four methods h must be less than 0.25. 
ШІ values of h < 0.25 , the Hamming method yields the 
 асест ассигасу, followed by the Fourth, Second, and 
ишта Order Adams'. The Milne method is a little better 
than the Second Order Adams'. Thus for accuracy the Hamming 
method ranks first, then the Fourth Order Adams, Milne, 
Euler, Second Order Adams, Third Order Adams, and lastly 
the Nystrom methods. 

For problems in the class of ODE it is recommended that 
h must be chosen less than 0.25 or smaller if valid numerical 
solution and accuracy are desired. Then use Hamming as the 
numerical method to solve the ODE. Again the Fourth Order 


Adams method is the second best choice. The choice of h can 
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TABLE 13 


TRUE SOPUTITONI VALUES rOR OpDE II 


Double Precision 


x Yexact 
1. КЕСЕНЕ 
ОЯ -4,931505:1071 
de 8.488724 101 
4. 1.410446 
5. 6.752620-10 + 
6. ESOS ROS 
72 -1.410888 
8. -8.438582:1071 
9. soon О 
10. 1.585092 





TABLE 13-A 


ОВЕ at x= 100 


ABSOLUTE ERROR РОВ Р-С-Г, РС ТГ, P C T l o p БЕКІ 


7590512" 


„5186358, 


1720949: 


70953907 


„2 ЛА 


‚160642 


131.604044 





-69510710 


.178410:10 


502092 710 


Double Precision 


‚634916°10 “17.950519:10 “|-1.764591 


.921730:10 9|5.455366:10 “| -6.947361 


-4 =] 


4.234564:10 =ОШ 


"515 214692 -100.834152 


-1123.407199 -353.009679 


‚482274 157.762746 -1081.991965 


4.69 5.34 221b 


116 





TABLE 13-B 
ODE at x = 10.0 


ABSOLUTE ERRORS FORMS Р-С-У ЕЕ 


Double Precision 


and РОС МЕРИ 


E HA 24 
5 |-1.20446-10 7 .760418: 4,448046-10 | -2.334506' 
Вт. 4837310” ‚192068. 3.610822:10 ?|-6.282905- 
32 | 5.531148: 10^ "|-2.331998: 2.961989°10 “| -1.347495- 
ДЕА 7. 53901 10 "У |-1.552622 БЕО СЕЕ тар 
И 686296- 10? |-7.673571” 2.016912 ‚600424. 
E 657253-10 2 |-6. 7065351: ILS S ens ‚268373. 
6. 2902132 44.871662 80.187459 ‚401846 

CPTS 5.26 ‚79 8.56 4.96 


ll; 








be formally stated as |hC| < 0.25, since the value of h to 


be used is directly dependent upon Ci- E Loy). 


ШІ ODE III 

Table 14 shows the exact solution values for ODE III from 
КЩ = 1.0 to x = 10.0. From the data of Table 14-A it is clear 
mat all four methods performed quite well to a good degree 
of accuracy for all balues of h. But for best accuracy and 
least computing time the Milne Method stands out followed by 
the Euler, Hermite, and Nystrom methods in that order. 

Erom Table 14=B, again. 1t 15 observed that any ot the last 
inr methods yield a valid numerical solution, to a certain 
Ewce of accuracy. For an excellent accuracy, however, the 
Hamming method is the best choice. Its range of accuracy is 


from 10 12 5 


up to 10) for values-of hīrangins М толу = 
to 1/h = 2. The Milne method compared to these last four 
methods ranks second only to the Hamming, though its com- 
puting time is a little smaller than Hamming's. Thus for 
excellent accuracy, the order of choice is an follows: 
Hamming, Milne, Fourth Order Adams, Third Order Adams, Her- 
mite, Euler, Second Order Adams, and lastly the Nystrom 


method in solving ODEs belonging to the class of ODE III. 


Ши ODE IV 

mue solution values are listed in ante el Set solace 
to x = 10.0. Actual data obtained in Table 15-A showed that 
all four P-C sets are good numerical methods for ODE IV. But 
1f a choice is to be made the Milne's accuracy is far greater 


than the other three methods with computing time only 0.16 


J Po 





TRUE SOLUTION VALUES FOR ODE III 


TABLE 14 


Double Precision 


119 


.414709* 
-0929745 
.411200: 
.568024* 
156890862 
.794154: 
3698609: 
EOS 
.121184* 


.440211* 


exact 


10 
107 
107 
10 
10. 
Oe 
ON 
10. 
107 


її 





TABLE 14-A 
ABSOLUTE ERROR FORE CI, PCIR C IIT An 


ODE III AT x=10.0 


Double Precision 


TI 


Id E .099:10. пао 50 


LI 


.845247* ПОО .814645: -288789 


' 70 3738: -2604110 .531438* .031052* 


.08984 5" OO > . 756454: .247845: 


. 365210: .596440*:10. .528314* .594418* 
.755574* -31713310 .856889* .260049: 
174787 .071599:1]0 .910558* .154229* 


2 АС 297 
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ПРБОЉОТЕ ERROR- FORSR ЕТ ТАРА А E ТА 


.01370'10 У 


.256927'10_ 
.970409: 
‚ 368642. 


1509286” 


TABLE 14-B 


ODE AT x = 


10% 


0 


Double Precision 


5.025175: 


790258: 


912798; 


.721024: 


.408691: 


. 767408 ` 


840913: 


1020 


. 375282: 


. 2870650: 


. 5919060. 


.588546: 


2850155” 


1040 


and P-C-VIII 


3.683:10 1l 


az 
121656 30 iO 
.649671:10. 
NA 7 On 


‚805833°10 


.360227:10 


„25 


4 








TABLE. ls 


TRUE SOLUTION VALUES FOR ODE IV 
= Double Precision 


X у 


ехаст 
ШЕШ) 3.894003 
200 5.659155 
ey 3.300505 
4.0 2.945055 
5.0 2576543 
6.0 22251509 
по ДАЗ 
8.0 1.624023 
oa №: 528155 
10.0 1.149189 
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TABLE 15-A 


ODE IV at x = 10.0 


ABSOLUTE ERROR FOR P C I, PC THT Ee Pi s o P G DD) 


Double Precision 


E 

‚544900 VL ‚598339: ‚575814: 
„839350. . .840047: . 826764: 
ee eos.) ‚339871° ‚278726 
926581: .97060*10 * |-2.872267- .889080- 
,176812- ,063982:1079|-1.119506:- 128099 
‚709013- ‚024899107 ‚260573- 418111: 
‚885932: ‚791009107 КЫ 668370: 
‚83 ‚67 6.74 
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sec longer than the next best choice which is Nystrom. After 
Nystrom, the Euler method is preferred; then ene serie 
Table 15-B data again demonstrates that any of the last 
four P-C sets can be used for solving ODE IV. Milne method 
compared with these methods competes with the Hamming in 
accuracy and has an edge in computer time (3.83 secs < 4.45 
secs). But this fraction of seconds difference is overcome 
Бу the demonstrated better accuracy of the Hamming in every 


a toh- 271. 


Step of the On from h = 2 Thus if 
preference is to be made, Hamming is the most likely choice 
ene best method for solving ODE IV class ot problems. 
The Milne, Fourth Order Adams, Third Order Adams, Nystron, 
Second Order Adams, Euler, and Hermite are the order of 


choices following the Hamming method. 


Б ОПЕ V 
Table 16 shows the true solution values for the range of 
integration considered. From the data of Tables 16-A and 


16-B, it is obvious that any method can be used to solve 


6 


ODE V and obtained accuracy up to minimum of 10 ^ and 


-14 


maximum of 10 Thus in comparing these methods, accuracy 


Exsteria must not be the greatest concern. By studying 
closely the behavior of the errors as h increases, it was 


noted that the Milne, Hamming, and Fourth Order Adams 


а 2:5 


exhibited decreasing errors from h = 2 to h = 2 for 


Milne, and from h = 220 to h = Да for both Hamming and 
Fourth Order Adams. Analyzing further the Hamming and 


Fourth Order Adams methods it could be seen that the Fourth 
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TABLE 15-В 
ODE IV at x = 10.0 


ABSOLUTE ERROR FORT CV, P-C-VI P CVII and c DD 
Double Precision 


1.28499:10 ? |-1.61:10 12 


025906: -2.690:10 11 


Ша 


то 282200: ‚174763. 3747-10 7° 


Ji 9 


NOS .884380* IS .20867:10 


9 7 


95941510 „131593 ames. 225495 210% 


E5192 107 .359162* ‚9200077 ° Ро 


E702 890 310. .625329* .006535* „201592510 


‚46 #02 2052 








TABLE 16 


TRUE SOLUTION VALUES FOR ODE V 


Double Precision 


x y 


exact 

О 1.020204 
2.0 1.040832 
50 J. 061515 
4.0 1.083472 
5.0 1a 105541 
6.0 109023152 
7.0 Jal51558 
SU 1.175759 
20 J 100555 
10.0 1.224744 
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TABLE 16-A 
ODE V at x = 10.0 


ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 
Double Precision 


ЕТ олоо оо | e o 
ТОО me w 
-609166* o s s mu 
.436058: EL Pornos ICT EET 
.960094* CEA o | Sun š 
978142: emo ое ІТ 
866820: 0408-107 .644850:1079|-1.478586:107” 
ЕЗ: Se 3.38 





14 


* X 
0 means error less than 10 


ШЕ 





ТАВ,Е 16-В 
ODESVoat xe) 


ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, 
Double Precision 


и 


and P(C YIII 


4.0°107 emote шо 1.6°107 

eons КОО ОООО | г а она E 

As e УЛО о [ооло RUP eens 
De o ота ооа 
Ое ao o e пао 
a ^ re ee ое 
,1744:10 19 | -7.649028:10 9 |-7.798880:10 9 |-2.7289:10 ^ 1? 

6.31 10.0 10.0 7.78 








Order Adams had its optimum value of-h at h < 27% while that 
of Hamming achieved its optimum h value at h < даа Thus 
the Hamming method has a higher value of h for its range of 
excellent accuracy. This property gives the Hamming method 
a slight edge over that of Milne, and the Fourth Order Adams. 
However if both accuracy and computing time are considered 
the Milne method is likely to be the choice. Therefore for 
the class of ODE V, the recommendation for the best method 
to be used will be most likely dependent upon the particular 
interest: whether accuracy and wider range of h values are 
the primary concern, or whether accuracy and least computing 
time is the criterion. For the former criteria the Hamming 
method serves the best purpose while for the latter the 
Milne method offers the best solution. However, as was 
noted earlier, if no criterion is involved but the interest 
is just to solve the problem the easiest way, the self- 
Ба пр Euler method is recommended. For realistic pur- 
poses however the following order of choice is highly 
recommended: Hamming, Milne, Fourth Order Adams, Third 


Order Adams, Hermite, Nystrom, Euler, and lastly Second 


Order Adams. 


Fee ODE VI 

True solution values are listed in Table 17 from x = 1.0 
to x = 10.0. Results from Table 17-A clearly showed that h 
must be chosen to be quite small to have a valid numerical 
solution. For the series of h values chosen none of the 


methods exhibited valid numerical Solution trom 1/h = 32 up 
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TABLE 17 
TRUE SOLUTION VALUES FOR ODE VI 


Double Precision 


< exact 
1.0 1.732050 
ди 2.236067 
3.0 2.645751 
4.0 3.0 

5.0 3.316624 
6.0 S OUS 
7.0 | 3.872983 
8.0 4.123105 
9.0 4.358898 
10.0 4.582575 





ПРАДУ ВАДЕ лана 


ODER Гас Хх ео 


S PSOLUTE ERROR FOR P< G-] РОСТОВЕ ТЕГ, ande C M 


Double Precision 









-50.844 MENOS -47.611 









20.0993 СИБИ - 0970061 





-214.085 6.986 133.800 


-436.640 21222700 





= 519.198 









23900292 





osos LED IGNIS 













219420. 155 210574 ES TES 







21095585 28979.47? -14737333 





205 ‚40 





БО? 





151 


- 32.764 


= оО 


18035158575 


592% 10% 


-5060.609 


“27792590 


-4096.870 


Zell 





То h = де The Milne method provided a valid numerical 


solution starting at h = 289 


and smaller, but the other 
methods still are unreliable. 
From Table 17-B, only the Hamming and Fourth Order Adams 


а апа 


method showed a valid numerical solution for h Z 
smaller, all others need much smaller h than the lowest 
value, 274 chosen. Again оп the average comparing the Milne, 
Hemming, and the Fourth Order Adams, for use in solving 


E the Hamming is the first 


ODE VI with values of h 2 
choice, Milne second and the Fourth Order Adams. As in 

the case of ODE II, the h values are governed by the maximum 
. magnitude of C = Е (х,у). If C is large then h must ibe 


enosen to be very small to satisfy the bounds for stability 


as had been established before. 


G. ODE VII 

True solution values are shown in Table 18 from x = 1.0 
to x = 10.0. Table 18-A showed that the Euler and Nystrom 
provide the accurate numerical solutions desired while the 
Milne and Hermite are unreliable. By analyzing the behavior 
of the roots of the Milne and Hermite methods shown in 
Table 18-B, it was observed that both exhibited unstable 
Seutions as evidenced by the uniform growth of the error 
in one direction whilethat of the Nystrom and Euler had 
decreasing error as the range of iteration increases.  Be- 
tween the Euler and Nystrom, the former produced mere еси 
rate results though it was a fraction of a second longer in 


computing time than the Nystrom. 
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TABLE 17-B 
ODE VI at x = 10.0 


I POOLUTE ERRƏDSFEFOBBRE C-V  P-G-VI sP C VI A 


Double Precision 


И 0 -9,371'10 7 


1.073'10" о бо В 


4.256 1 012695 A 5357 


4.556 -547.874 2457839 - 354: 261 


= 202.070 55 063 198.05] -100. 57] 


БОТ 2015: 7854 2322.89 ZO 


800.00] - 521022 3621.438 7807 7305 


6.68 216 9 
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TRUE 


105 


TABLE 18 


SOLUTION VALUES FOR ODE VIT 


Double Precision 


0 | 
9 
0 5) 
9 
0 9 
9 
9 
9 
9 
0 D 


154 


1999097 


Y exact 


.615941: 
040275 
.950547: 


OZ 


„999877 
. 999983. 
. 999997. 
. 999999. 


„999999 


1071 


1071 


107? 


1071 


= | 


Ae 


1071 


1071 


по 





TABLE 18-A 
ОБЕ УТ ato т! 


ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 
Double Precision 


175.005 72 
7091970 
ASUS ; . 0655106: 


‚5577-1019 


„000472: „4540,10 ‚511858: 


10 


8764-10 322982: ‚345210 10 


о 


.58311:10 "| 3.314797 .46305:10 ^? .187207* 


.07753:10^?| 1.429417 .593088*10^ .852296* 


4.42 4.75 4.14 „03 
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TABLE 18-B 
ODE at h = 0.5 


ABSOLUTE ERROR БОР Р-СІІ,РСІ NIE NE 


Double Precision 


So 


you? 


. 540" 


7200 


554; 


MOS 


ENS 


8527: 








Results from Table 18-C showed that the Пас КОШЕ оке 
produced valid numerical solutions with Hamming and Fourth 


Order Adams competing for best accuracy. As h increases 


zu 25 


from h = 2 to h = 2 


-4 


the Hamming method is most precise 


=] 


but from h = 2 to + 2 the Fourth Order Adams method 


showed better accuracy. Computing time for the fourth Order 
Adams method is less than the Hamming method. The Euler 
method compared with these two methods as h increases, 


"к D 


=] 


started at lesser accuracy from h = 2 but 


and yielded 
zu 


maintained its accurate results up to h = 2 
a much more accurate numerical solution at h = 2 It also 
had the least computing time. Thus for classes of ODE 
belonging to ODE VII, the recommended order of choice of 
methods is as follows: Euler, Nystrom, Fourth Order Adams, 
Hamming, Third Order Adams, and lastly the Second Order 
Adams. The Milne and Hermite are not considered and should 


Moi De used for this particular type of ODE since they are 


both unstable. 


PRODE VIII 

lable 19 presents the true solution values for ODE VIII 
from х - 1.0 to x = 10.0. Data analysis from Tables 19-A 
ШЕШІ |І9-В indicates that any of the eight P-C sets can be 
msed to solve ODE VIII if no specific criterion for accuracy 
is needed, since each method yields. good and valid numerical 
solutions. However, for purposes of comparison, the Hamming 
method clearly stands out to be the best in terms of accuracy 


with Milne coming in second. The Milne method provides the 


253% 





ТАВЬЕ 18-С 


(DELE Ex SEU О 


ESSOLUIE ЕКО РОВ ВО, РУ, РОЗИ, 


Double Precision 


п A a eee 


7.665862 ° 


ER 574510- 


(052571: 


1.049942: 


SES 


Ооо 


ADOS 


561650: 


ЛОИ 


,008471: 


‚3506166 


mA 


AU 


5. 0668010 О" 


ZOGEN 


D I EOS 


Порто аи 


119715 10 
29290099910. 


.033756 "10. 


a20 


and P-C VITI 


0 


-10714 


о 
-71598510 
7270 13 10) 
.570674°10 > 


‚82 








TABLE 19 


MAME: SOLUTION УАТПЕ5 FOR ODE VITI 


Double Precision 


X Y exact 
1.0 Е 

2 DAT 

5.0 D 

4.0 4.691641'10 1 
5.0 zo Se PORE 
6.0 В ОО а 
7.0 | — 1.928970 

8.0 2.689507 

9.0 1.510013 
10.0 5.804096:10 1 
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57125277 


552090" 


1050551 ° 


. 505628" 


8521252” 


520272: 


1051157” 


TABLE 19-A 


ODE at x 


10. 


0 


Double Precision 


.17509:10. 
.498490 10. 
.440709:10 


1450298: 10. 


.793708 10 


КОПТ БЛ Ол 


‚63 


2) 


8 


7 


6 


5 


5 


140 





611170: 


204757: 


8506620 


.254094* 


. 557320: 


„1152458 


270577075 


К 


ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, 


and P-C-VI 


Lug 


«5958457 


"03954 


.499001` 


58050808 


641152. 


2705 5 


‚866206 


6.48 





TABLE 19 B 
ODE at x = 10.0 
ABSOLUTE ERROR FOR P-C:V,; P-C-VI BS VITARA El 


Double Precision 


10 


9,902758:10. 8983 107 


3.303524:10 9|1.886624:10^ .53331:10 ^? 


7.843814* .498027- 11559975089 


| zi 


.006735-10. -1.359301: ЕЗ 5 ТӨМ” ‚563035°10 


6 E 


2112250 * 105 UNE AA: 1249289: -197783 10 


Кос то ет .783513- O 


4 = 


"95617110 SG DAE „783513 „11317010 


"o'l IT 5. 09 








least computing time with Hamming second. Thus for problems 
of the type ODE VIII, the Hamming method is the most aia 
recommended, followed by the Milne, Fourth Order Adams, 
Euler, Hermite, Nystrom, Third Order Adams, and lastly the 


Second Order Adams. 


ime ODE IX 

True solution values are listed in Table 20 for ODE IX 
ШЕШ (пе range of integration x = 1.0 to x = 10.0. From 
Tables 20-A and 20-B it was observed by analyzing the 
meomlts that all eight methods provided valid numerical solu- 
tions. In terms of better accuracy on the average and the 
least computing time the Milne method is the best candidate 
Сибел Палта пе coming in next, followed by the Fourth Order 
Adams, Hermite, Third Order Adams, Euler, Nystrom and 
lastly the Second Order Adams. Thus if choice is to be 
made for solving classes of ODE belonging to ODE IX, the 
Milne method is highly recommended with Hamming as the next 


alternate. 


МОРЕ Х 

ls differential equation is an example of a controlled 
variable problem, where one variable is expressed as a func- 
tion of x and substituted into the differential equation of 
the Other variable to form a single first order ODE. The 
Original two variable equations are: 


у! = 1.38y, - 0.81y, do S 


Да до мод аве, (2-75) 


142 





TABLE 20 


PRUE SOLUTION VALUES. FOR, ODE ITX 
Double Precision 


К Yexact 
j. O 2.000555 
20 2.249705 
5.0 4.179509 
4.0 9.462527 
pea 10.633343 
6.0 17.564095 
O 42.421552 
8.0 50.806495 
9.0 74.608406 

10.0 186. 463649 
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ТАВЪЕ 20-А 
ODE IX atx = и 


ABSOLUTE EPROR FOR p G J РС p cI P n ес 


Double Precision 


amm 5 








ТАВ,Е 20-В 
ODE IK at x = 10.90 


ABSOLUTE ERROR FOR P-C-V, P-C-VI, P-C-VII, and P-C-VIII 
Double Precision 








with initial conditions 
y (0) = -2.9995 (2-76) 
Y, (0) = 4.0010 (2-77) 


and the analytic solutions given by 


0.0005 ox - зе 0: 3Х 


N 


У} (2-78) 


0.001 е 2Х + 4e 9: 5X 


|| 


У (2-79) 


Snte the analytical solutions of both differential equations 
are known, it is possible to treat one variable as a function 
of x and substitute it to the differential equation of the 
gener variable to form a single first order ODE which now can 
ge olved by the P-C sets. Choosing y, as the known function 
ШЕП Dy (2-79) and substituting it in (2-74), the resulting 
equation 1s ODE X with initial condition given by (2-76) and 
eier solution by (2-78). Note that у; remains a function 

of the single independent variable x and as such changes as 

x changes. The other variable Y 1 forms a first order dif- 
ferential equation in the standard form У] = f(x,y). The 
pent that rs being illustrated here is that this concept can 
be applied to problems like rate of chemical reaction, falling 
Beatles, aircraft or ballistic missiles flight wherein time 
variable is the most important consideration. All other 
Epuables can be"expressed as Tünctronsgosmerhe sqm MN 
pendent variable time and given constants or initial condi- 
tions. In such cases the problem can be reduced to a single 


differential equation and the methods considered herein can 
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be applied (time (t) variable аз х). This concept does not 
discount the fact that these P-C methods can be extended to 
solve simultaneous differential equations with some modifi- 
cations in the algorithms. However, since this extension 
is not relevant to the purpose of this particular paper in 
analyzing and comparing the predictor=corrector methods аз 
applied to a variety of ODEs, it is not considered here, 
but will be mentioned as a further field of study. 
Returning now to the solution of ODE X, Table 21 shows 
the true solution values from x = 1.0 to x= 10.0. From 
Table 21-A the Milne seemed to produce the best accuracy 
We by analyzing further the behavior of these errors in 
step-by-step integration, which is shown in lable -ICE Rt 
is noted that the errors of the Milne method are increasing 
mm a uniform fashion such that as the range of integration 
is increased the error grows while the true solution values 
as shown in Table 21 decrease. Thus eventually the error 
will overcome the solution. The same is true for the 
Hermite method, while for the Euler and the Nystrom the 
C Or decreases as the range of integration increases thus 
providing a valid numerical solution. The behavior of the 
Milne and Hermite methods is expected since for ODE X it is 
easy to see that C < 0, and as such the Milne and Hermite 
methods are unstable. Between the Euler and the Nystron, 
the Euler method is an easy choice both in accuracy and 


computing time. 
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TABLE 21 


ТКОЕ SOLUTION VALUES FOR ODE X 


Double Precision 


x Yexact 
ies -2.222429 
2. К ООЛ; 
3, -1.219708 
n 
4. -9.035826'10 
E 
5. -6.693904°10 
й 
o» -4.958966:10 
7. ОО 
n 
8. КОК ОЛЕ, 
ОП 
9. ТО поделена 
10. -1.493612 10? 
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TABEE 212 
ODE atx = 1029 


ABSOLUTE ERROR FOR P- C-I, PCI PCI., nik €n 


Double Precision 


HE 
‚839775. L 52: -9,818242:1074 
‚887212. ,3459:107 584152: 3.913399. 10, 
.893782+ .606138* 832304: -1.558160:10 2 
‚149738 - .286345* 793736: разсада елша 
757046: Е „903244 SS 
800289 Ear 095989: -1.288607 
784363" 436505: ‚ 350331 - -11.095378 
o 10.0 9.83 
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JABLE 21 D 


ODE XK. with he 0s 


ABSOLUTE ERROR FOR P-C-I, P-C-II, P-C-III, and P-C-IV 
Double Precision 
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From Table 21-C, all the last four methods exhibited 
stable numerical solutions. For accuracy the Hamming method 
offers the best choice, with Fourth Order Adams coming 
next, though the computing time for Hamming is a little 
longer than the Fourth Order Adams, = therefore tor probléns 
Ша (Пе class of ODE X, the order of recommended preference 
15: Hamming, Fourth Order Adams, Third Order Adams, Euler, 
Nystrom, and the Second Order Adams as the last choice. The 
Milne and Hermite methods should not be used for this par- 


ticular type of problem. 


K. SUMMARY 


7 фо ћ = 273 and тапре 


For series of h values from h = 2 
e tegration up to x= 10.0, Tables 22, 22-A, and 22-8 
list the summary of results for the eight predietor Correeton 
sets considered, each run on the test ODEs I to X. From 
en se comparative results, it is clearly established that 
the best numerical methods in order of decreasing efficiency 
in general are: 

1. Hamming Modified P-C Set 

2. Fourth Order Adams -Mouliton P-€ Set 


Milne P-C Set 


> WW 


Third Order Adams-Moulton P-C Set 
eee E Шер 

6. Nystrom P-C Set 

7. Second Order Adams-Moulton P Gg en 


8. Hermite P-C Set. 
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ТАББЕ А G 
ODE X at x = ]0 0 
ABSOLUTE ERROR FOR P-G-V, P-C-VI, 


P C VII P Pp ае IET 


Double Precision 


E ЕЛ 
‚476336 о 
. 538981: 949594: 3801011 
,773099:1079|-5.207967: 912304: 064 ы 
.065215-107 .272440* 180575: K 
1708110” 630456: Бат, ,553614:10 ^ 
ШОО J o КОТ .628554* .293084:10 Š 
9.45 ‚98 5.45 
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TABLE 22 


SUMMARY OF RESULTS EORZODEZT ZTOZOBE ZEN 


то ће2 3 AT x-10.0 


FOR h=2 ' 








— 






Computing 

































ODE Stable P-O Sets in Average Ter Unstable 
Number Order of Preference Accuracy rd P-C Sets 
I 1. Hamming 1. Milne 
2. Fourth Order Adams 2. Hermite 
3. Third Order Adams 
EME 
5. Nystrom 
6. Second Order Adams 
11 1. Hamming l Hermite 
2. Fourth Order Adams 
3. Milne 
ИЕ CT 
5. Second Order Adams 
6. Third Order Adans 
7. Nystrom 
ШТ 1. Hamming Хопе 
2. Milne 
3. Fourth Order Adams 
4. Third Order Adams 
5. Hermite 
6. Euler 
7. Second Order Adams 
8. Nystrom 
ТУ 1. Hamming Да Мопе 
2. Milne 3. 
5. Fourth Order Adams 4. 
4. Third Order Adams Te 
ЕТОМ 3u 
6. Second Order Adams Dom 6.46 
7. Euler TUNE 4.39 
8. Hermite 1074 6.74 
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TABLE 22 
SUMMARY OF RESULTS FOR ODE V TO ODE VIII 
FOR h=27 


O A ио 








Computing 
Time 
in sec. 


Unstable 
P-C Sets 


s ES Ss = <-> zx was 


Hamming 

Milne 

Fourth Order Adams 
Third Order Adams 
Hermite 

hystrom 


Euler 


ОО -і О ол FP U N HH 


Second Order Adams 








. Hamming і : Еси 
Fourth Order Adams I : . Nystrom 
. Milne : . Hermite 


Third Order 
Adams 







. Second Orde 
Adams 















. Milne 


Hermite 


Ји Шет 
Nystrom 
Fourth Order Adams 
Hamming 
. Third Order Adams 


Second Order Adams 










CN Qn + (4 ҒӘ FF 


Hamming 

Milne 

Fourth Order Adams 
Euler 

Hermite 

Nystrom 

Third Order Adams 


Second Order Adams 


u о (n co co ve pP MII DN A A + 


Con СО YN + (А NH 
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ODE 
Number 


IX 


TABLE 22-B 
SUMMARY OF RESULTS FOR ODE IX AND ODE X 
БОҚ һ-277 ТО h=2"1 AT x=10.0 














Computing| Unstable 
Tıme 
in secs. 






Stable P-C Sets in 
Order of Preference 


Average 
Accuracy 





Milne 

Hamming 

Fourth Order Adams 
Hermite 

Third Order Adams 
Fil 

Nystrom 


СОО Nu CN «л A U N e 
-1 —1] DW со CN Cn Cn сл 


Second Order Adams 






l1. Milne 


2. Hermite 





Hamming 
Fourth Order Adams 
Third Order Adams 


Euler 









Nystrom 


С л > (4 н 


Second Order Adams 





The most important criteria applied iS accuracy mola 
computing time. It is obvious that the Hamming methoaidid 
not produce the least computing time as compared steer 
other methods especially with that of Milne and the Fourth 
Order Adams methods. This difference of less than a second 
in computing time is attributed to the fact that the Hamming 
method used additional computation for the truncation error 
to modify the predicted and corrected values for every 
iteration. Analyzing the computing time for each method 
“thin the range of their stability, the list below shows 


the order of increasing computing time: 


Average 
DENS CIE Computing Time 
jue AO LS CCS 
2. Hamming 5.60 secs. 
3. NyS trem 5.057°scec s; 
4. Fourth Order Adams 6.11 secs. 
S. ЕПТ 6.19 secs. 
6. Hermite 6.58 secs. 
7. Second Order Adams 8.40 secs. 
8. Third Order Adams 8.69 secs. 


ШИЕ сопіттпесй the idea that convergence ог а соттестст сес 
not necessarily assure the convergence towards the true 
solution values, y (x,) > but only to some definite value, 

ul This is best illustrated in the case of Hermite and 
Third Order Adams methods with 6.58 secs and 8.69 secs average 
computing time, respectively. By recalling from the order 


of efficiency of the P-C sets, the Third Order Adams ғалыс 
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fourth in accuracy while Hermite ranks eighth пп о ли око и 
to the fact that though the Hermite method converee Smich 
faster than the third order Adams, its accuracy is much 
poorer than the Third Order Adams. Clearly then it can be 
seen that the Hermite method is converging only to some 
definite value, E but not to the true solution, у(х), 
otherwise it should have much better accuracy than the Third 
Order Adams method, which took much longer to converge to 


afaetinite value y in each iteration. This analysis is 


пе е 
meaningful due to the fact that a standard mode of applica- 
tion is used, that of Bee. It should be noted however 
that rapid convergence is essential to accuracy as shown 

by the Hamming, Milne, and Fourth Order Adams methods which 


all rank well above the other methods in both accuracy and 


memnmputing time in general. 
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А, CONCLUSTONG 


The following conclusions can be drawn and recommenda- 
tions made from the analysis and numerical results obtained 


Mi this paper. 


1. Numerical instability results from extraneous solutions 
of the difference equations which bear no connection to 
the exact solution. The conditions of asymptotic (strong) 
and absolute stability can be seen clearly by reference 
tO the extraneous solutions. If in the limit h+- 0, the 
extraneous solutions vanish as n > o, then the method is 
Eemunoly stable and convergent. If, for values of h less 
than some ho» the extraneous solutions vanish as n >œ, the 
method is absolutely stable. Relative stability is a sig- 
wM acant concept for ODE where f > 0. The condition pro- 
vides that extraneous solutions will not grow more rapidly 
c Cay more slowly than the true solution. Thus, before 
a method is used, the characteristic roots of its related 


difference equation must be analyzed. 


2. The finite real negative stability bounds for the P-C 
meenods considered, as determined by numerical experimenies, 


are in reasonable agreement with theoretical predictions. 
As such experimental bounds could be used as a guide to the 


proper selection of a method to solve a particular оа оше 


3. The stability of a P-C set depends on both the predictor 


and corrector equations. Thus, when two P-C sets have the 
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same corrector but different predictors, choose the one with 
a better predictor. Likewise, when two P G Бес СКЕ 
same predictor with different correctors, choose the one 


with a better corrector. 


EER C sets with higher order truncation errors ao 
necessarily better. In choosing a P-C set,-the order of 


truncation error should not be the sole criterion. 


5. The convergence of the corrector formula will ensure 
that the sequence of approximations will converge to some 
definite value but not necessarily to the true solution 
values. As such, the fast computation time of a method 
oce not necessarily imply greater accuracy for the result- 


ing numerical solution. 


6. If the integration will involve a large number of steps, 


a stable method should be used. 


feet the function evaluation is lengthy, the range of inte- 
Ма оп large, and better accuracy and lesser cost of com- 
puting time are prime considerations, then the following P-C 
methods are recommended, based on the overall efficiency 
they have demonstrated in numerical experiments on a wide 
variety of different test ODEs: 

a) Hamming P-C Sets 

b) Fourth Order Adams P-C Set 

c) Milne P-C Set | 

а) Third Order Adams P-C Set 

e) Euler P-C Set 


INS ROME Бет 
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g) .Second Order Adams P-C Set 


h) Hermite P-C Set 


The listing indicates order of preference. However caution 
must be exercised not to use Milne and Hermite methods if 


а < 0. 
У 


8. For maximum accuracy to be achieved, select a step size 
h, Гог which the truncation error and the roundoff error have 


equal orders of magnitude. 


9. If high accuracy results with less roundoff error are 
Екей, (Пе ргосеайите should use double precision, which 
involves only slightly more computing time than single pre- 


cision in an IBM 360/67 computer. 
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XI. EXTENSIONS 


The following topics, which are direct extensions of 
this paper, are interesting areas of further research and 


study in the wide field of numerical analysis: 


ime Stability versus Accuracy 

A necessary condition for numerical integration is 
stability; that is,a stable value of h must be employed 
when a fnite stability boundary exists. Nevertheless, 
there is no guarantee that all values of h from zero to the 
limiting vlue will yield accurate .résults. Thus the 
question that must be faced is: Of what value is a method 
that is stable in a region where it is inaccurate? The 
answer to this question will show that the relationship 
between stability and accuracy is fundamental to the choice 


of h for a particular method. 


ШИ! Іпстеавсіпс the Stability Bounds 

It has been shown that in some methods the real stability 
bound is somewhat constrained. Two possible ways where the 
stability bounds can be increased are: 

a) By the use of aweighting factor which involves experi- 
menting with different values of the free parameters that 
are used in the method by sundetemmined) coci.1c lento, cea 
to Hamming's corrector formula derivation. 

b) By averaging, which means computing the new value of 


т after a finite number of steps (say fifty) of a particular 
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method, as the average of тие оу and a value called 


п+ 1 
Yari generated by another method. This should be done 


periodically. 


КО Predictor Correcto Netnods Interaction 

The interaction of the different methods in solving a 
Mr ticular problem involves use of the P-C sets as sub- 
routines, wherein one method will be used to provide the 
solution for smaller values of h, then another method will 
bemused to solve for larger values of h. This should be 
quite interesting since it has been shown that some methods 
exhibit better accuracy for small values of h, then become 
inaccurate for large values of h, while other methods main- 
pun their accuracy up to large values of h but might be 
prohibitive to use for small values of h due to their longer 


computing time. 


4. Solution of Systems of Differential Equations 

Revision of the algorithms to enable the solution of 
©% ems of differential equations and higher order differ- 
ential equations. This extension is straightforward. It 
involves only extra computational steps using the same 
formulas. Complexity arises from the need to use vectors 


and arrays for temporary storage. 


5. Practical Applications 
Applications to real-life problems such as heat flow 
problems, simple electrical circuits, force problems, rate 


of bacterial growth, rate of decomposition of radioactive 


FG? 


material, crystallization rate of a chemical compound ено 


of population growth, and so forth. 


0. Adjustment of the Step Size During the P C o Miron 

This involves monitoring the modified truncation error 
developed by Hamming. If the value is too large, then 
the step size is too large, and the calculation should be 
repeated with a smaller value of h, say h/2. Note that it 
winot necessary to go back to the beginning of the calcula- 
tion but only to the point at which the truncation error 
becomes too large. If the truncation error is too small, 
computation time is being wasted and h should be increased, 


say to 2h. 
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FLOW CHART FOR ALGORITHM 1 


INPUT TERMNNAL BOUNDARY, 
NUMBER OF ITERATIONS, 


CONVERGENCE TERM, NUMBER 
OF STARTING POINTS 





READ STEP SIZE AND 
CONTROL PRINTING 


ШБЕСІНІ INITIAL CONDITION 


END Yes 
OF INPUT DATA 
SET? 


No 


PRINT INPUT DATA SET 
CALL FOURTH ORDER 


RUNGE-KUTTA METHOD 










I5 
METHOD SELF- 
STARTING: 


TO GENERATE 
NEEDED ADDITIONAL 
SIARTING POINIS 
COMPUTE FUNCTIONAL 


EVALUATIONS OF 
STARTINGTPOTNIS 
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COMPUTE- PREDTCTED 
SOLUTION VALUE USING 
PREDICTOR FORMULA 


COMPUTE FUNCTIONAL 
EVALUATION OF PREDICTED 
SOLUTION VALUE 





COMPUTE CORRECTED 
SOLUTION VALUE USING 
CORRECTOR FORMULA 


COMPUTE MAGNITUDE 
OFIDIEFEFEBENCETSBETI EEN 


PREDICTED VALUE AND 
CORRECTED VADPUE 









СОМРАКЕ 
DIFFERENCE WITH 
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NUMBER 
OF ITERATIONS 
SATISFIED: 


SET PREDICTED 
SOLUTION VALUE 
EQUAL TO CORRECTED 
SOLUTION VALUE 





COMPU TEARS 
SOLUTION VALUE 


COMPUTE NUMERTCAL 
ERROR IN CALCULATION 


PRINT 
RESULTS 


UPDATE REQUIRED 
PARAMETER VALUES 


ADVANCE INTEGRATION 
POINT DY ATOTEP 
SIZE INCREMENT 












INTEGRATION 
BOUNDARY 
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FLOW CHART FOR ALGORITHM 2 





INPUT TERMINAL BOUNDARY, 
NUMBER OF TTERATIONS; 
CONVERGENCE TERM, NUMBER 


OF STARTING POINTS 
READ S Ter SIZE 
AND CONTROL PRINTING 


SPECIFY INITIAL CONDITION 


END 
OF INPUT DATA 
SEE, 














Yes 






PRINT INPUT DATA SET 






CALL FOURTH ORDER 











RUNGE-KUTTA METHOD 





IS 


AENA ЕШТЕ 
STARTING? 










TO GENERATE 






NEEDED ADDITIONAL 








STARTING POINTS 






COMPUTE FUNCTIONAL 
EVALUATIONS OF 
STARTING OLN TS 
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COMDUTE PREDICIED 
SOLUTION VALUE USING 
PREDICTOR FORMULA 


COMPUTE MODIFIED 
PREDICTED SOLUTTON 
VALUE 


COMPUTE FUNCTIONAL 
EVALUATION OF MODIFIED 
PREDICTED VALUE 


COMPUTE CORRECTED SOLUTION 
VALUE USING CORRECTOR 
FORMULA THEN COMPUTE MODI- 
FIED CORRECTED SOLUTION VALUE 


COMPUTE MAGNITUDE OF 
DIFEERENCETBETWEEN MODI- 
ЕВЕ ЮР ОЕ ШЕЕ И ЕЕЕ 
MODIFIED CORRECTED VALUE 
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COMPUTE TRUE 
SOLUTION VALUE 


COMPUTE NUMERIAAL 
ERROR IN CALCULATION 


PRINT 
RESULTS 


UPDATE REQUIRED 


PARAMETER VALUES 





ADVANCE INTEGRATION 
POINT BY A STEP SIE 
INCREMENT 


ERMINA 
INTEGRATION 
BOUNDARY 
EXCEEDEDZ 


No 


о 
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КОМВЕК 
ОВ ВЕЕРОМ 5 
ШЕШЕ Ш/2 















SET MODIFIED 
PREDICTED VALUE 
EQUAL TO MODIFIED 
CORRECTED VALUE 






Y'--Y 


ODE: 


Jt Jt 3 3E 3E 3t 36 3€ 36 36 36 3E 3E 3t 2E 3€ 2E 3C 3€ 3€ 3C 36 ЗЕ ЗЕ ЗЕ Зе ЗЕ ЗЕ 36 


TO SOLVE ANOTHER ODE SIMPLY CHANGE 


PROGRAM ONE 
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EULER PREDICTOR-CORREC TOR METHOD. 
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INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 
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o~x~n<x= 
MODO Il 
"ИН И d e 
<-<-п-<» 


TEST FOR CONVERGENCE 

LECDELY=EPS? 305303525 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 
al es 20,20,30 


ou 


co TO 11 

COMPUTE TRUE SOLUTION VALUE 

YEXACT=EXACT(XO) 

COMPUTE ERROR IN CALCULATION 

ЕККОК=ҮЕХАСТ-ҮС 

TESI IF DESIREO POINT CF PRINTING TS REACHED 


IF((N/ICON*ICONJ).EQ.N) ИАІТЕ (6, 1002) ХО, ҮР, ҮС, ҮЕХАСТ, 
ТЕККОК 


UPDATE REQUIRED PARAMETER VALUES 
YO=YC 
FP=FC 
40 CONTINUE 
LOOP BACK TO READ NEXT. SET OF INPUT DATA 


UJ 
© 


ЛЛ SETTUSEDZIYZTFGSRK TER =", 
15 OF EULER PREDICTCR- CORRECTOR 
Dr ҮП 1 р ҮС 


' ЕККОК !//12X; о 44X y 
ХаҒ15./7,2Х,ЕҒ15.7:2Х,Ғ15.7,2Х,Е15.7,2Х, 


O am P a pm 

„< ММТ 

O Mix> не 
X "Dn. NO 
PNA se 

О О ~ 
a Dim 

O в Соч 


— 
< 
== 


< 


sN IPN 
——— QM 


М РСТ (ХМ. УМ) 


mann 
EME 


ON ЕХАСТ(ХХ) 
EXP (XX) 


mm” 
ЕС 
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*INPUT DATA SET USED* 
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H = 0.5000000 
ICON - С 
**RESULTS OF EULER PREDICTOR=COPRRECIO Ев 5 
Х УР ҮС ҮЕХАСТ ЕВКОК 
0.0000000 1.0000000 
1.0000000 0.3588867 0.3634033 0.3678794 0. 0044 761 
2. 0000000 О. 503212 0.1325720 051352555 0.0027633 
3.0000000 0.04 776 13 0.0483635 0.0497871 0.0014236 
4. 0000050 0.0174238 0.0176435 0.0183156 0.0006722 
5.0000000 0.0063564 0.0064365 0.0067379 0.000301% 
6.0000000 0.0023189 0.0023481 0.0024788 0.0001307 
1.0000000 0.0008459 0.0008566 0.0009119 0.0000553 
8.0000000 0.0003086 0.0003125 0.0003355 0.0000230 
9 . 0000090 0. 0991126 0.0001140 0.0001234 0.0000094 
10.0000000 0.0000411 0.0000416 0.0000454 0.0000038 
*INPUT DATA SET USED: 
H - 1.0000000 
ICON = 1 
F RESULTS OF FULERTPREDICTOR-CORRECTOR METHOD 
X YP YC YEXACT ERROR 
0. 0000000 1 .0090000 
1.0000000 0.2500000 0.3750000 0.3678794 -0.0071206 
2. 0000000 0.15625 00 0.1 718750 0.1253353 -0.0365391 
3.0000000 0.0507813 0.0683594 0.0497871 -0.0185723 
4.0000000 0.0253789 0.0300293 0920153156 =—0¿011713/ 
5.0000000 0. 0095825 0.0122986 0.0067379 -0.0055606 
6.0000000 0.0044327 0.0052910 0.0024788 -0.0028122 
1.0000000 0.0017519 0.0021987 0.0009119 -0.0012868 
8.0000000 0.0007731 0.0009362 0.0003355 -0.0006007 
9.0000000 0.0003156 0.0003919 - 0.0001234 -0.0002685 
10.0000000 0.0001361 0.0001660 0.0090454 -0.0001206 





PROGRAM TWO 


: үз=-ү 


ООЕ 


MILNE PREDICTOR-CORRECTOR METHOD., 


TO SOLVE ANOTHER ODE SIMPLY CHANGE 


INITIAL CONDITIONS ANO FUNCTION SUBROUTINES. 


Y(O)=1. 


WITH 


3 


f 
О 
THE FOLLOWING 


el ll ==80.<6. 
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«(OX LOFFOO-OOO-OrLa2 Daek 
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GZ >ü Oé OZu.,— u Does CO9ED 
WO TY SOM DRO Wh K ОД и 5 О 3: 
ROS Ak OF LLM ODU p- Ėja he je je pe Us 
IS DO IE) 3t 
LU Iu uc uiii E mulu a Se 
EEZIEOITISTZZTSSDEZE Lim ie 
fe fe ee DR PE ОР Е LL U. U о 2: +: 
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u u сооз 





Y ОКШ << >ш ж >— > 


а. со 


EXACT E 


RK UA. 
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jt 3E 3€ 3E JE 3E 36 3€ 3E 36 3C 3C 36 36 JC 36 OE 3C XC C 3E 3E 3E 3 3€ 3E 303€ 303€ 


о ооооеооомосоооооооооооооомоооооооомоеоо 


INITIALIZE INPUT CONSTANTS 


AX=10.0 
X=2 
еи 4220001 


ХМ 
МА 
ЕР 
K= 


READ STEP SIZE AND CONTROL PRINTING 


КЕАО (5,100) 


Н, ІСОМ 


5 


SPECIFY INITIAL CONDITION 


OOO Фо 


ОО 

° О 
O > 
нии 


>< > > 


TEST FOR END OF INPUT DATA 


IF(ICON-0) 4514537 


PRINT INPUT DATA SET 


OOO фо 


NEXT SOLUTION POINT 


COMPUTE 


OW 
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ооо ООС 


ооо ооо обо 


AAD 


AAN 
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бок 
Ou 


UJ 
O 


40 


МС< (ХМАХ+Н/2. ) ИН 

Х1<ХО+Н 

X2=X0r2.*H 

X3=X04+3 . *H 

Х=Х 0+4. *Н 

CALL RKUTTA(XO+YR,Y19Y2,Y3,K,H) 
FI=FCT(X1,Y1) 

Ғ2-ЕСТ(Х2,Ү2) 

F3=FCT(X3,Y3) 

по 0 N=“ , NC 
YP=Y04+4.*H/3.*(2.*F1-F2+2.*F3) 
ЕР-ЕСТ(Х,ҮР) 

YC=Y 2+H/3.*(F2+4.*F3+FP) 
DELY-ABS(YC-YP) 


TEST FOR CONVERGENCE 

IEEBELY-EPSI 30730715 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 
ІҒ(М-МАХ) 20,20,30 

FP=FCT (X, YC) 

МЕМ+ 1 

GO TO 11 

COMPUTE TRUE SOLUTION VALUE 

ҮЕХАСТ-ЕХАСТ(Х) 

COMPUTE ERROR IN CALCULATION 

ERROR=YEXACT-YU 

IES? DESTIRED POINT ОР PRINTING 15 REACHED 


A 00 XYP YC УЕХАСТ, 


UPDATE REQUIRED PARAMETER VALUES 


F1=F2 
F2=F3 
F3=FP 
YO=Y1 
Y1=Y2 
Y2=Y3 
Y3=YC 
X=X+H 
CONTINUE 
LOOP BACK TO READ NEXT SET OF INPUT DATA 
СО ще 5 
FORMAT (F10.7+110) 
1T(////41X+* “INPUT DATA SET USED*=*//43X+'H =*, 
О. 1/ ^3X, * ICON =*,110) 
AT[//29X, t XSRESULTS. OF MILNE PREDICTOR- -CORRECTOR 
METHOD*-*'//]10Xx.' X 135 X4, Y ! YC 
7 ҮЕХАСТ DR ERROR 4//12Х 10.7, 44Х, 
` ee ҒЕ15.1,2Х3Ғ1І6. ТЭ2ХТЕГО ТҰ2ХУКТО УТ ОАО 
ЕМО 
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FUNCTION FCT(XN,YN) 


FCT=-YN 
RETURN 


END 


TION EXACTEXK) 
Та 


МС 
АС 
TU 
D 
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*INPUT DATA SET USED? 


= 0.5000000 


2 


PREDICTOR-CORRECTOR METHOD** 


YC 


0.1353100 
0.0496315 
0.0181091 
0.0064749 
0.0021262 
0.0004233 


-0.0003513 
-0.0013280 


YEXACT 
1.0000000 
0.60653 06 
0.367894 
0.2231301 
О 203558585 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


INUNDA TAS EII: 


H 
ICON 
**RESULTS OF MILNE 
X YP 
0.0000000 
0.5000000 0.6067709 
1.0000000 0.3681710 
1.5000000 0.2233955 
2.0000000 0.1385594 
3.0000000 0.0509259 
4. 0000000 0.0183159 
5.0000000 0.0061856 
6.0000000 0.0015049 
Т. 0000000 -0.0005340 
8.0000000 -0.0017356 
9.0000000 -0.0028179 
10.0000000 -0.00%1233 
H 
ICON 
RESULTS OF MILNE 
X YP 
0.0000000 
1.0000000 0.3750000 
2.0000000 0.1406250 
3.0000000 0.0527344 
4.0000000 0.0468752 
5.0000000 0.0147570 
6.0000000 0.0102883 
7.0000000 0.0014462 
8.0000000 -0.0005506 
9.0000000 -0.0039288 
10.0000000 -0.0064589 


= 1.0000000 


1 


ЕККОК 


-0.0002403 
~0.0002916 
-0.0002654 


0.0000253 
0.0001555 
0.0002066 
0.0002630 
0.0003525 
0.0004886 
0.0006868 
0.0009704 
0.0013734 


DREBITLTUÜR-CORRECTERZPRETHODFZ 


ҮС 


0.0164931 
0.0051922 
0.0002441 
0.0005433 


-0.0009222 


0. 0010040 


—0. 0017257 
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YEXACT 
l «0000000 
0.36 18794 
02713592593 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 


-0.0071206 
-0.0052897 
-0.00294 73 


0.0018226 
0.0015457 
0.0022346 
0.0003686 
0.0012576 


-0.0008806 


0.0017711 
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ОРЕ: 
FUNCTION SUBROUTINES. 


PROGRAM THREE 
TO SOLVE ANOTHER ODE SIMPEY CHANGE 


NYSTROM PREDTETOR-CORKECTIOR METHOD: 
Ү(0)-1. 
INITIAL CONDITIONS AND 
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ООО Ооо 


ооо ооо ооо 


eere 


OOO 
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o 
© 


1002 


1004 


МС+ ( ХМАХ+Н /2. ) /Н 
X1=X0+H 

X=X0+2.*H 

CALL РКОТТА(ХО, ҮР, Ү1.Н) 
ЕІ= РСТ (Х1,Ү1) 

eno N=2 + NU 
YPzYO-*2.*H*Fl 
ЕР=ЕСТ(Х,ҮР) 

YC=Y 1+H/2.*(F1+FP) 
DELY=ABS(YC-YP ) 


TEST FOR CONVERGENCE 

Е СРЕСУ= ЕРОЛ ЗО: Зо Б 

TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 
1Е (М-МАХ) 2020, 30 

FP=FCT (X,YC) 

M=M+ 1 

GO TO 11 

COMPUTE TRUE SOLUTION VALUE 

УЕХАСТ+ ЕХАСТ (Х) 

COMPUTE ERROR IN CALCULATION 

ERROR=YEXACT-YC 

TEST 1F DESIRED POINT OF PRINTING IS REACHED 


NO WRITE(6,1002) X,YP,YC,YEXAGT, 


1 

UPDATE REQUIRED PARAMETER VALUES 

YO=Y1 

Yl-YC 

Fl-FP 

ХЕХ+Н 

CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

СО ТО 5 

STOP 

FORMAT (F10.7,110) 

FORMAT (////41X+** INPUT DATA SET USED**//43X, 'H =", 
1F10.71/43X+=* ICUN =*,110)] | 

FORMAT (//29X_**RESULTS OF NYSTROM PREDICTOR-CORRECTOR 

IMETHOD = © /7/ LOK! Х 4,5Х,! ҮР Ay, ШЕ i 
ПЕР ЕЕЕ ЭХ." ЕККОК “//12Х,Ғ10.7,44Х, 

ЕОҚМАТ(%04%,9Х,Ғ15.7,2Х,Е15.7,2Х,Е15.7,2Х,Ғ15.7,2Х; 
1Ғ15.7) 

ЕМО 

УОВЕОЈТТМЕ ККОТТАСХХ , У У, У], НН) 

СКІзННжЕСТ(ХХ,ҮҮ) 

CK22HH*FCT(XX*HH/2.0,YY*CK1/2. 0) 

CK3-HH*FCT OXX*HH/2.0, YY*CK 2/230) 

CK4zHH*FCTOXXTHH,YYtTCK3 

Yl-YYt(CK1*2.0*CK2*2.0*CK3t*CK41/6.0 

ХХЕХА+НН 

УЕХАСТ+ЕХАСТ (ХХ) 

ERROR=YEXACT-Y 1 

НАІТЕ(6, 1004) ХХ, У], УЕХАСТ , ЕВЕСВ 

Е АТОМНОЕ БОНА О s ала 


NCTION ЕСТ( ХМ, УМ) 


ЕМО 


ОМ ЕХАСТ(ХК) 


МСТ 
ХАСТ-ЕХР(-ХК) 
Е К 


Le Wy OF tL 


МСТ] 
А = 
TURN 
D 
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*INPUT DATA SET USED* 


*RESULTS OF NYSTROM PREDICTOR-CORRECTOR METHOD * 


X 
0.0000000 
0.5000000 
t. 0000000 
2. 0000000 
3. 0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


*RESULTS OF 


X 
0.0000000 
1.0090090 
2.0000000 
3.0000009 
4.0000000 
5.0090000 
6.0000000 
1.0000000 
8.0000000 
9.0000000 
10.0000000 


H = 0.5000000 
ICON = 
YP YC 

0.6067709 
0.3932291 0.3636068 
0.1444499 0.1298206 
0.0516075 0.0463004 
0.0184066 0.0165119 
0.0065643 0.0058886 
0.0023410 0.0021000 
0.0008349 0.0007489 
0.0002977 0.0002671 
0.0001062 0.0000952 
0.0000379 0.0000340 


УЕХАСТ 
1.0000000 
0.6065306 
0.3678794 
021555353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
070005255 
0.0001234 
0.0000454 


*INPUT DATA SET USED* 


H = 1.0000000 
ICON = 
NYSTROM 
YP YC 
0.3750000 
0.2509000 0.1093750 
0. 0625000 0.0156250 
0.0468750 -0.0053594 
-0. 0015125 020078125 
0. 0097656 -0.0041504 
—-OVOOSTUST 2202002177753 
0.0046387 -0.0005798 
-0.0045166 -0.0003052 
040028381 0. 00005 72 
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YEXACT 
1.0000000 
0.3678794 
0421553355 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.9001234 
0.00004 54 


ЕККОК 


-0.0002403 


0.0042726 
0.0055147 
0.0034867 
0.0018037 
0.0008494 
0.0003 788 
0.0001630 
0.0000684 
0.0000282 
0.0000114 


PREDICTOR-CORRECTOR METHOD ж 


ERROR 


=0. 007/1206 


0. 0259603 
0.0341621 
0.0241750 
0.0145504 
0.0066291 
0. 0031091 
0.0009153 
0.0004286 


—0.0000115 
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СОЕ A 


TO SOLVE ANOTHER ODE SIMPLY CHANGE 


INITIAL CONDITIONS AND FUNCTICN SUBROUTINES. 


PROGRAM FOUR 


HERMITE PREDICTOR-CORRECTIOR METHODS 
Y(0)=1. 


WIYH 
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= = 


10.0 
0.000001 


ICON... 
YE XACT 
ERROR 
ЕСТ 
EXBET__ — 
RKUTTA 
AAA AA AAA MR A z £ 
AX= 
Х=2 
READ STEP SIZE AND CONTROL PRINTING 
READ(5,100) H, ICON 


5 


INITIALIZE INPUT CONSTANTS 


SPECIFY INITIAL CONDITION 


ХМ 
МА 
ЕР 


3€ 36 30 3€ 36 3€ 36 26 3€ 3€ 3E 3C 3E 3€ 36 30 36 36 3C E 3E CIE dE E CE 303 


5 


QOOQUOQOUOQOQQOOOQOQOQOQUQOOOQOUOOQQOQOQOOQGOOOOQO0 OOO OOO 


ИЕ 


TEST FOR END DF INPUT DATA 
CUMPUTE NEXT SOLUTION POINT 


IF(ICON-0) 45,45,7 
PRINT INPUT DATA SET 





ооо MAM 


Ооо обо ООО ООС 


AAN 


H/2.)/H 


Ни Ни 


AA s: 
) 


YltH*(4.*F1*2.*FO0) 


ii F044. *F1+FP) 


O <NZONTOXx>xZz 
mio o 9 II OF2 O D> IM 


r Un U 


TEST FOR CONVERGENCE 

IF(DELY-EPS) 30,30,15 

TEST IF MAXIMUM NUMBER QF ITERATIONS IS SATISFIED 
ТЕ (М-МАХ) 20,20,30 

ЕР+ЕСТ (Х, УС) 

M=M+1 

GO TO 11 

COMPUTE TRUE SOLUTION VALUE 

ҮЕХАСТЕЕХАСТ(Я) 

COMPUTE ERROR IN CALCULATION 

ERROR=YEXACT-YC 

TEST IF DESIRED POINT CF PRINTING 15 REACHED 


DU LEONE *ICON).EQ.N) WRITE(6,1002) X,YP,YC, YEXACT; 


De 
Ou 


о 
© 


UPDATE REQUIRED PARAMETER VALUES 
Ү0-Ү1 
NE 
FO=F1 
Fl-FP 
X=X+H 
40 CONTINUE 
LOOP BACK TO READ NEXT SET OF INPUT DATA 
GO TO 5 
45 STOP 
100 FORMAT(F10.7,110) 
1000 FORMAT(////41X,'*INPUT DATA SET USED*'//43X;, "Н =", 
1F10.7/43X,* ICON -*, 110) 
1001 FORMAT(7/29X, * RESULTS OF HERMITE PREDICTOR-CORRECTOR 
1METHOD *1//10X, Хх 1,5%, x С ; 
15. УЕХАСТ 1,5%, "ERROR "//12А;Ё10.7,44&Х, 
1002 ЕОВМАТ(*'0*,9Х,Е15.7,2Х,Е15.7,2Х,Е15.7,2Х,Р15.7,2Х,‚ 
1F15.7) 
END 
SUBROUTINE RKUTTA(XX,YY Yl, HH) 
CK1=HH*FCT(XX, YY) 
CK2=HH*FCT(XX+HH/2.0, YY+CK1/2.0) 
CK3-HEECT XXe HH/2 20 2 YY «CR 2/210) 
CKA-HH*ECT(UXX*HH , YY*CK3) 
Y1=YY+(CK1+2.0*CK2+2.0*CK3+CK4)/6.0 
ХХЕХХ+НН 
YEXACT=EXACT(XX) 
ERROR=YEXACT-Y1 
WRITE(6,1004) XX,Y1, YEXACT , ERROR 
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FUNCTION FCT(XN, YN) 
ЕСІ--ҮМ 
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Х 
0.0000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6. 0000000 
1. 0000000 
8.0000000 
9.0000000 

10.0000000 


RESULTS OF 


X 
0. 0000000 
1.0000000 
2.0000000 
3. 0000000 
4+. 0000000 
5, 0000000 
6. 0000000 
7. 0000000 
8. 0000000 
9. 0000000 


*INPUT DATA SET USED* 


ERROR 


-0.0002493 
0.0002818 
0.0003963 
0.0006143 
0.0009343 
0.9014098 
0.0020879 
0090031057 
0.0046165 
0.0068610 
0.0101963 


METHOD * 


ERROR 


~0.0071206 
0.0057054 
-0.0234295 
0.0275871 
0202217095 
0.0863936 
-0.1470038 
0.2470207 
-0.4164615 


H = 0.5000000 
ICON = 2 
*RESULTS OF HERMITE PREDIC TOR-CORRECTGOR METHOD = 
YP YC YE XACT 
1.0000000 
0.6067709 0.6065306 
063593160 0.3615916 0.3618794 
0.1296880 0.1349390 0:1353853 
0.0435692 0.0491728 0.0497871 
0.0106042 0.0173813 0.0183156 
-0.0043826 0.0053371 0.0067379 
=0.01390225 0.0003908 0.0024788 
=. 023417 9972 2029921378 0.0099119 
-0.0358057 -0.0042811 0.0003355 
=0. 0535 12217 = 0927 00636 0.0001234 
-0.0797584 -0.0101509 0.0000454 
<INPUT DATA SET USED* 

H = 1.0000000 

ІСОМ - 1 
HERMITE PREDICTOR~CORRECTOR 
ҮР ҮС ҮЕХАСТ 
1.0000000 
0.3750000 0.3678794 
0.0090000 ОРО 00 0313533953 
01620950 0.0732166 0.04978 71 
—-0.2105594 099092115 0.0183156 
0.3834665 0.0599074 0.0067379 
0.6344536 -0.0839149 0.0024788 
1.0731830 0.1479157 _ 0.0009119 
-1.8065000 -0.2466852 0.0003355 
340441850 0.4165850 0.0001234 
=5,.1255590 591014316 0.0000454 


10.0000300 
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HAMMING PREDICTOR-CORRECTOR METHOD. 


TO SOLVE ANOTHER ODE SIMPLY CHANGE 


-1. 


NIYH Y(0) 


FUNCTIGN SUBROUTINES. 


INITIAL CONDITIONS AND 
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INITIALIZE INPUT CONSTANTS 


XMA X 


10.0 
0.000001 


2 


READ STEP SIZE AND CONTROL PRINTING 


КЕАО (5,100) 


Н, ТСОМ 


SPECIFY INITIAL CONDITION 


а 
OOO OVO 


TEST FOR END OF INPUT DATA 


IF(ICON-O) 45,45,T 


PRINT INPUT DATA SET 


OOO YOO 
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COMPUTE NEXT SOLUTION POINT 
NO nU ы 


0,ҮК,Ү1.Ү2.ҮЗ,К,Н) 


>> -< -< -< 7> 
w w ">Ç 


AUN eme 
р 
IND O—-M OUpM—-— 


<< we ee = 
= Ó) 


Ө ОО Ae ЕЗ) 
ТЕСЕ ЕТ 


.X*Y3-Y1)43.*H/8.*( FP*2, F3-F2) 
*ET 
) 


11 


r" H < Se 
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We 


TEST FOR CONVERGENCE 
IF (DELY-EPS) 30,30,15 
TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 
15 ТЕ (М-МАХ) 20,20,30 
20 ЕЕ. FCT (X,YC) 
=M+] 
GO TO 11 
COMPUTE TRUE SOLUTION VALUE 
УЕХАСТ + ЕХАСТ (Х) 
COMPUTE ERROR IN CALCULATION 
ERROR=YEXACT=YC | 
TEST IF DESIRED POINT OF PRINTING IS REACHED 


1F((N/TICON=ICON).EQ.N) ЯКІТЕ(6,1002) Х-ҮР,ҮС,ҮЕХАСТ, 
1ЕККОК 


UPDATE REQUIRED PARAMETER VALUES 


Fl=F2 
F2-F3 
F3-FP 
YO=Y1 
үшү? 
Y2=Y3 
ҮЗ=ҮС 


X=X+H 
40 CONTINUE 
LOOP BACK TO READ NEXT SET OF INPUT DATA 


е С 


= егп || Jt 
се г 

x rd 

> 


Ма. SET USED*"//43X,'H = 
F HAMMING PREDICTOR~CORRECTOR 
,52Х,! ҮР 19X ҮС ' 

‚> ERROR 1128, Flo. 1144Xy 

,Ғ15.7,2Х,Ғ15.7,2Х,Ғ15.7,2Х,Ғ15.7,2Хх, 


VO 
ec 
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ION EXACT(XK) 
EXP (XK) 


М 


СТІОМ ҒСТ(ХМ,ҮМ) 


FUNCT 
EXACT 
RETUR 
END 





*INPUT DATA SET USED* 


*RESULTS OF HAMMING PREDICTOR-CORRECTOR METHOD = 


Х 
0.0000000 
0.5000000 
1.0000000 
1.5000000 
2.0000000 
3.0000000 
4 ,0000000 
5. 0000000 
6. 0000000 
7. 0000000 
8. 0000000 
9.0000000 

10.0000000 


Х 
0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6 , 0000000 
7. 0000000 
8. 0000000 
9.0000000 


Н = 0.5000000 
ІСОМ - 2 
ҮР ҮС 
О „606 7 709 
0. 3681710 
022253855 
0.1385594 0.1355409 
0. 0494538 0.04992 89 
0.0184383 0.0183953 
0.0068499 0.0067794 
0.0025701 0.0024990 
0.0002748 0.0009215 
0, 0003766 0.0003401 
0.0001496 0.0001256 
0.0000619 0.0000465 


ҮЕХАСТ 
1.0000000 
0.6065306 
0.3618194 
0.2231301 
0% 19592952 
0.0497871 
0.0183156 
0.0067379 
0.002.788 
0.0009119 
050002355 
0.0001234 
0.0000454 


“INPUT DATA SET USED* 


H = 1.0000000 

ICON = 1 
ЖЫБИШЕТ5 OF HAMMING PREDICTOR-CCRRECTOR 
YP YC YEXACT 
1.0000000 
9.3 150000 0.3678794 
0.1406250 0.1353353 
0.0527344 0.04978 71 
0.0468752 0 0190868 0.0183156 
=O. 0199183 0.0056581 0.0067379 
0.0245465 0.0057350 0.0024788 
-0.0516121 —0.0008928 0.0009119 
0.0793315 0.0053606 000059355 
—0.118#7?798 -0.0069360 0.0001234 
0.1838089 0.0110968 0.00004 54 


10.0000000 
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ERROR 


-00.0092403 
-0.0002916 
-0.0002654 
-0.0002056 
-0.0001418 
-0.0000806 
-0.0000415 
—0. 09000202 
-0.0000096 
-0.0000046 
0200090072 
-0.0000011 


METHOD * 


ERROR 


-0.0071206 
—0.0052897 
-0.0029#+73 
=0. 0007112 

0.0010798 
-0.0032592 

O. 0018047 
-0.0050251 

0.0070594 
-0.0110514 
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SECOND ORDER ADAMS PREDICTOR=CORRECTOR R MET DDD 


TO SOLVE ANOTHER CDE CHANGE 


CONDITIONS AND FUNCTION SUBROUTINES. 
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TEST FOR CONVERGENCE 
IFLEDELVY-EPS17305539,15 
TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 


ТЕ (М-МАХ) 20,20, 30 
ЕР+ЕСТ (Х, УС) 


GO TU 11 

COMPUTE TRUE SOLUTION VALUE 

ҮЕХАСТ-ЕХАСТ(Х) 

COMPUTE ERROR IN CALCULATION 

ERROR=YEXACT-YC 

TES! IF DESIRED POINT OF PRINTING IS REACHED 


IF((N/ICONZICON).EQ.N) WRITE(6,1002) Х,ҮР,ҮС, ҮЕХАСТ, 
1ERROR 


ALA SET USED**//43X,'H =>, 
S OF ADAM2 PREDICTOR-CORRECTOR 
АЕ ҮР 15X3! YC '. 
3 ERROR “"//12Х,Ғ10.1,44Х, 


ке 
m n ДД Ce a at) 


:YY4Y1,HH) 


Ү%СК1/2.0) 
ее 
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FUNCTION ЕСТ (ХМ, УМ) 


FCT=-YN 
RETURN 


RETURN 
END 


END 


ТОМ ЕХАСТ(ХК) 
a 


FUNCT 
EXACT 
RETUR 
END 


КӨҢ 





Х 
0.0000900 
0.5000000 
1.0000000 
2. 0000000 
3.0000000 
4. 0000000 
5. 0000000 
6.0000000 
7. 0000000 
8. 0000000 
9.0000000 


ERROR 


-0.0002403 
0.0044048 
0.0105122 
0.0148087 
0.0171469 
0.0182921 
0.0188207 
040190556 
09121571 
0.0192002 


*INPUT DATA SET USED* 
Т = 0.5000000 
ICON = 2 
**RESULTS OF ADAM2 PREDICTGR-CORRECTOR METHOD** 
ҮР ҮС УЕХАСТ 
1.0000000 
0.6067709 0.6065306 
0.4016927 0.3634746 0.3678794 
0.2968013 0.1248231 0.1353353 
0.2556222 0.0349784 0.049787і 
0.2401217 0.0011687 0.0183156 
0.2342887 -0.0115542 0.0067379 
0.2320936 -0.0163420 0.0024788 
0.2312676 -0.0181437 0.0009119 
0.2309567 -0.0188217 0.0003355 
0.2308398 -0.0190768  0.0001234 
0.2307958 -0.0191728 0.0000454 


10.0000000 


Х 
0. 2000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


*INPUT DATA SET USED* 


H 


ICON 
**RESULTS OF ADAM2 


YP 


0.3750000 
0.3125000 
0.343 7500 
0.3281250 

459909315 
90292202513 
0.3339844 
0.3330078 
0.2334961 
0.3332520 


1.0000090 


1 


0.0192182 


PREDICTOR~CORRECTOR METHOD** 


YC 


0. 10056025 
—0. 0312500 
-0. 0996094 
moe CIN 5 
2031499023 
= Је Гагра! 
-0.1624 756 
-0.1645508 
=), 16550789 
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YEXACT 
1.0000000 
0.3678794 
О 13553535 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0009454 


ERROR 


= 02 0077205 
0.03377238 
0.0810370 
0.1179250 
0.139550 
021523239 
0% 15397750 
0.1628110. 
0.1646742 
0.1656643 
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TEST FOR CONVERGENCE 
ТЕ (ОЕГУ-ЕР5) 30,30,15 
TEST ТЕ MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 
15 ТЕ (М-МАХ) 29,20,30 
20 ЕР-ҒСТ (Х,ҮС) 
=M+1 
GO TO 11 
COMPUTE TRUE SOLUTION VALUE 
ҮЕХАСТ-ЕХАСТ(Х) 
COMPUTE ERROR IN CALCULATION 
ERROR=YEXACT-YC 
TEST IF DESIRED POINT OF PRINTING IS REACHED 


ЕР IE EN МЕІТЕ(6,1002) Х,ҮР, ҮС, ҮЕХАСТ, 


UJ) 
e 


DATE REQUIRED PARAMETER VALUES 


<<< ттт 
TORK Te ~ 


n uut 


40 CONT INUE 
LOOP BACK TO READ NEXT SET OF INPUT DATA 


ја A 

O OMO 
о оо 
Fe 


pauses SET USED**//43X,'H = 
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po 
po 


0%,9Х,Ғ15.7,2Х,Ғ15.7,16Х,Ғ15.7,2Х,Ғ15.7) 


Cea 


о 
© 


FUNCTION FCT(XN,YN) 


CTION EXACT( XK) 
EXP (-XK) 


105 





*INPUT DATA SET USED* 


H 
I CON 


"RESULTS OF ADAM3 


X 
0. 2000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000900 
6. 0000000 
7. 0000000 
8. 0000000 
9.0000000 
10.0000020 


MI? 


0.60677 09 
0.3681710 
0.2630881 
0.0949250 
020532505 
0.001118 
0.0014410 
0.0005048 
0.0001768 
0. 0000619 


- 0.5000000 


2 


PREDICTOR-CORRECTOR METHOD** 


ҮС 


0.1307259 
0.0457278 
0.0160214 
0.0019661 
0.0006887 
0.0002412 
0.0000845 
0.0000296 


YE XACT 
1 -0000000 
0.6065306 
0.3678794 


‚0.1353353 


0.0497871 
0.0183156 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


*INPUT DATA SET USED* 
= 1.0000000 


H 
ICON 


**RESULTS OF ADAM3 


Х 
0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
1. 0009000 
8. 0000000 
9.0000000 
10.0000000 


ҮР 


0.2750000 
0.1406250 
0.1888022 
0.0217022 
0.0785821 
-0.0525024 
-0.0479046 
0. 05171729 
0.0527026 
-0.0540226 


1 


ЕККОК 


-0.0002403 
-0.0002916 


0.0046093 
0. 0040593 
0.0022942 
0.0005126 
0.0002232 
0.0000942 
0.0000389 
0.0000158 


PREDICTOR-CORRECTORTMETHOD S% 


YC 


00.0454788 
0.0021872 


-0.0024490 
-0.0057783 


0.0015025 


—0. 0918528 


0.0029584 


-0. 0020290 
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ҮЕХАСТ 
1. 0000000 
0.3678794 
05.19 53859 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 


-0. 0071206 
—0. 0052897 


O. 0043082 
00.0161284 
O «0091869 
O. 0082570 


-0. 0005906 


0.0021883 


—0.0028350 


0.0020744 
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TO SOLVE ANOTHER ODE CHANGE 


PROGRAM EIGHT 
1. 
INITIAL CONDITIONS AND FUNCTION SUBROUTINES. 
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Ү(0) 


Y'!'=-Y WITH 


ee 


+ —! + 
x < + 
„ш ы зе 
ка. ос — + 
Or > LLJ ~ ш ш 3€ 
+ | I ч (> ul (C 257 t 
X QU u. €: =. PUES = Sc + 
+ SOOM WY XI TU) = > “ 
та ха >? Ze oL > [EX Ти! Us $ 
шо OX Xx O еч Or ш kk + 
*xOXNX>O RRR ш < OU Ur | Al + 
жано Lu гіт Z м ТЕТЕ з 
+ ныт ч < <1 о м-н Бор а ш + 
жол о ou Сір-ггіш Ш = нах 2 > Z + 
иг IET Zuim UO> О ш d 
От OZrer Fr С 2-45» 9 «x + 
3% ON ee VO == 2 Uu O E ша Ww 3t 
3 1 ош O uù u ONA ra КЕ 00 45 
+ у> > ШО O „Dr CSOD B uu 
ЈЕ се С AO | <\1ООСО<— 22% 
ко фт < UJ  ашуоцэе>г---ш-. Ш О 
+ DZ Qua: III SS = x 2 2 Wat 
90 мых zu J + ШНУРА e оь 
¿Our D Ot «i Or-RZONXO>TF rr X 
3 me Ощш8> > >р-2--о-и Dv (956 
го MOD Fr Rw -»OruwuFOQO о ошшек 
Ооо гого WW ара 2 пресни 
Ее Оғы Ju Welle Sel 2 ОШ Осо | ср ређе ЈЕ 
j€ 2uJu—O. <(ОО>З-ОР-ООК"-"ОООРеОК-СС0о DIe xk 
j O T 2:00 Оов сб. ш о NIF? 
LA JUHO Se Шш к eœ u. Ob 3$ 
Зе = е тате Аш DOLOS EE 
скен ООСС ОСЕ ШД Ооо ООоОО с + 
+ О ч LU CO «1 £ — th OO сны DO 
+U ZIO ноа. ни оао әш <Б ББ ош) 
х OW} =, О. со DOVOD OFX 
36 = (О бу = UW U UI U oad U oad U UU UU U EE Z LLULLU JU Z Z Z —JUJLL)3E 
tacu ттт: 
кш со PrreNrerrroZzee0mnu uu mt 23 
$ Owu «y шиш | | | | 1. | | Жы 
ЈЕ (О на > | | | || | | | | | | | х 
Xo OC ew | | | + 
FO >= | | % 
% О»>--<9 | | < % 
ноша | | | | | | | Oe | — ы +z 
3 << C) <( lÍ < Z >< | | || = < O tO e % 
Е о T го xo j < о +» 
YX O г ха «ча. © lv ше O XxX Y + 
ос О ОМ < 026 Ц 2: > = ШІ € хш ци шах 
х оош» гш x 
Я ч 
jt 3€ 3€ 3€ 3€ X6 3€ 36 30 3€ 3€ 36 3€ 3€ 3€ 4€ 26 3€ 3€ 3€ XE 3€ 2€ 9€ 30 306 3€ 3€ 3€ X 
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10.0 
0.000001 


2 
READ STEP SIZE AND CONTROL PRINTING 


INITIALIZE INPUT CONSTANTS 
READ(5,100) H4ICON 


XMAX 
SPECIFY INITIAL СОМОТТТОМ 


5 


QOO омео 


оо 

e eO 
Q> 
ини 


>< >~ >- 


Ш 


TEST FOR END OF INPUT DATA 
COMPUTE NEXT SOLUTION POINT 


ТЕ (ТСОМ-0) 45,45,7 
PRINT INPUT DATA SET 
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ТЕЗІ EGR CONVERGENCE 
LELBELY-ER5) 30,30,15 
TEST IF MAXIMUM NUMBER OF ITERATIONS IS SATISFIED 


15 IF(M-MAX) 20,20, 30 
20 с 
= М+ 


CaL CI] 


COMPUTE TRUE SOLUTION VALUE 
30  YEXACTsEXACT (X) 

COMPUTE ERROR IN CALCULATION 

ERROR=YEXACT-YC 

TEST IF DESIRED POINT OF PRINTING IS REACHED 

IF((CN/ICON*S ICON) .EQ. N) WHRITE(G6,1002) X,YP,YCS, YEXACT, 

1 ERROR 

UPDATE REQUIRED PARAMETER VALUES 

YO=Y1 

Yl-Y2 

Y2=Y3 

Y3=YC 

FO=F1 

Fl=F2 

F2-F3 

F3-FP 

X=X+H 
40 CONTINUE 

LOOP BACK TO READ NEXT SET OF INPUT DATA 

С@ [0 5 
45 STOP 
100 FORMAT(E10.7+110) 
1000 FORMAT (////41X,**INPUT DATA SET USED*'//43X,'H =", 

ТЕТ 43%." TOON ТО) 
1001 FORMAT (//29X, ***RESULTS OF ADAM4 PREDICTOR-CORRECTOR 

lIMETHODS*'//10X,* X ',5X,' ',5Х%,' Y m 

1547. YEXACT — ',5X;,' ERROR  "//L2X ELO АН С 
1002. FDRRAT (*0* ,9X1F15- T 2X 4 F15 T: 2X4 F15, 7) 2X5 F 15-1 2X 

END 
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ООТТМЕ ККОТТА(ХХ,ҮҮ,Ү1,Ү2,ҮЗ,КК,НН) 


17 


*CK3+CK4)/6.0 


HH 


ХХ, ҮҮ, ҮЕХАСТ,ЕККОК 


ЕХАСТ(ХХ) 
УЕХАСТ-УУ 


+2. ОЖСКЗ+СКА ) /6. 0 
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1004 ҒОКМАТ(%0%;9Х,Ғ15.7,2Х,Ғ15.7,16Х,Ғ15.7,2Х,Ғ15.7) 
NCTION ЕСТ (ХМ, УМ) 
199 
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“INPUT DATA SET USED* 


200 


H = 0.5000000 
[СОМ = 2 
*#RESULTS OF ADAM4 PREDICTOR-CORRECTOR METHOD** 
X YP YC YEXACT ERROR 
0.0000000 1.0000000 
0.5000000 0.6067709 0.6065306 -0.0002403 
1.0000000 0.3681710 0.3678794 -0.0002916 
1. 5000000 0.2233955 0.2231301 -0.0002054 
2.0000000 0.1397457 0.1352788 0.353353 0.0000565 
3.0000000 0.0512678 0.0495746 0.0497871 0.0002124 
4. 0000000 0.0187947 0.0181669 0.0183156 9.0001 488 
5.0000000 0.0068879 0.0066569 0.0067379 0.0000810 
6.0000000 0.0025232 0.0024393 0.0024788 0.0005394 
7.0000000 0.0009245 0.0008938 0.0009119 0.0000180 
8.0000000 0.0003387 0.0003275 0.0003355 0.0000079 
9.0000000 0.0001241 0.000 12.00 0.0001234 0. 0000034 
10.0000000 0.0000455 0.0000440 0.0000454 0.0000014 
*INPUT DATA SET USED* 
H = 1.0000000 
ICON = 1 
"RESULTS OF ADAM4 PREDICTOR=CORRECTOR METHOD = 
X ҮР ҮС ҮЕХАСТ ЕККОК 
0.0000000 1.0000000 
1.0000000 0.3 150000 0.3678794 -0.0071206 
2.0000000 0.1406250 0513933537 79 00922 71 
3.0000000 0.0527344 0.0497871 ~0.0029473 
4. 0000000 0.0744628 0.0149522 0.0183156 0.0033634 
5.0000000 0.0091045 -0.0007950 0.0067379 0.0075320 
6.0000000 0.0319238 -0.0004662 0.0024788 0.0029450 
7. 0000009 -0. 0306429 -0.0027268 0.0009119 0.0036387 
8.0000000 0.0368799 9.0015700 0.0003355 —0.0012345 
9. 0000000 ~0.0442679 -0.0027738 0.0001234 0.0028972 
10.0000000 0.0550162 0.0028112 0.0000454 -0.002 1658 
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